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APPLICATION  OF  NON-ORTHOGONAL  CURVILINEAR 
COORDINATES  TO  CALCULATE  THE  FLOW 
IN  TURBOMACHINES* 


Chen  Nai-xing** 

(Marine  System  Engineering  Research  Institute) 


In  the  first  part  of  the  paper,  the  formulas  of  calculating  gradients  of  a  scalar, 
divergents  and  vorticities  of  a  vector  in  non-ortliogonal  curvilinear  coordinates  are  pre¬ 
sented  by  veetor  analysis.  With  aid  of  these  relations  we  then  obtained  basic  aero- 
thermodynamic  equations  governing  relative  steady  flow  of  a  nonviscous  fluid  in  a  tur¬ 
bomachine.  General  non-orthogonal  coordinates  are  suggested  for  solving  0-equation 
of  three-dimensional  turbomachine  flow.  The  potential  equation  for  three-dimensional 
flow  calculation  is  obtained. 

IP-equations  (Wu’s  equations),  expressed  in  term  of  general  non-orthogonal  coords 
nates,  of  two  kinds  of  stream  surfaces  are  discussed  in  the  third  part. 

In  the  last  part  of  this  paper,  three  forms  of  velocity  gradient  formulas  are  pre¬ 
sented. 


1.  NON-ORTHOGONAL  CURVILINEAR  COORDINATE  SYSTEM 


The  Publication  of  Comrade  Wu  Zhonghua's  paper  [2]  has  attracted 
attention  from  everywhere.  By  using  the  tensor  method,  he  derived  for 
the  first  time  the  fundamental  equations  for  3-dimensional  flow  In 
turbomachines  in  non-orthogonal  curvilinear  coordinates,  which  stimu¬ 
lated  the  development  of  a  design  method  for  turbomachine  aerodynamics. 
We  started  our  work  in  1973.  The  publication  of  reference  [2]  has  also 
stimulated  and  inspired  our  work. 


* 

This  paper  has  been  presented  at  the  Second  National  Engineering 
Thermophysics  Conference  in  Hangzhon  in  November,  1978. 

*§ 

Currently  at  the  Mechanics  Institute,  Academic  Sinica. 
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X* 


Figure  1.  Non-orthogonal  Figure  2.  Non-orthogonal 

curvilinear  coordinates.  curvilinear  micro-element. 

Key:  1 — constant. 

In  a  non-orthogonal  curvilinear  coordinate  system,  the  position 

12  3 

of  a  point  M  may  be  represented  by  the  3  coordinate  values  x  ,  x  ,  x  . 
The  3  unit  vectors  un ,  u~,  u,  are  taken  in  the  tangential  directions 

-L  ^  j  12  3 

along  the  increment  directions  of  the  3  variables  x  ,  x  ,  xJ  at  the 
12  3 

point  M.  u  , . u  ,  uJ  represent  respectively  the  3  unit  normal  vectors 
perpendicular  to  the  3  surfaces  of  the  non-orthogonal  curvilinear 
hexahedron  with  volume  W 

1  p  3 

The  partial  derivatives  of  the  radial  vector  rQ(x  ,  x  ,  xJ) 
may  be  written  as 

drjdx'  ™  HiUit  (f-  1,2,3).  (1.1) 

where  is  the  length  of  drjdx'  ,  called  the  Lame  coefficient. 

They  may  be  obtained  from  the  following  equation: 

w?  “  (drjdx1)1  —  (dx/dx‘Y  +  (dy/dx'y  +  (dx/dx')*.  (1.2)  /111 

Any  vector  C  decomposed  along  the  2  sets  of  3  non-coplanar  vectors 

12  3 

u-^,  u2,  and  u  ,  u  ,  uJ  may  be  written  as 

C  -  C'a,  +  C*o,  4-  C*o„ 

C  —  CtB*  +  CjU1  +  CjB*.  (1.3) 
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(1.4) 


The  relationship  between  and  C1  is 

C,  —  [C'cosCa,,  Oi)  ■+■  Clcos(aj4  o2)  4-  C’co s(a„  a01/«»(aS  a,) 


As  is  well  known,  a  vector  pointing  in  the  direction  of  the 
fastest  change  in  a  scalar  *  with  a  magnitude  equal  to  the  derivative 
in  this  direction  is  called  the  gradient  of  $  .  The  derivative  along 
an  arbitrary  direction  m  is  the  projection  of  the  gradient  (v$  )  of 
*  in  this  direction.  According  to  this  definition,  it  is  not 
difficult  to  derive  the  formula  for  the  gradient  as: 


— 


1 


d<t> 


cos(a1,U|)  H\dxl 


a1  4- 


d<p 


a2  4* 


d<P 


cot  (a2,  Bj)  Htdxl  cos  (a*,  a,) 


(1.5) 


The  flux  of  a  vector  C  per  unit  volume  through  the  surface  of  an 
infinltisimal  volume  around  a  point  is  called  the  divergence  of  the 
vector,  i.e. 


V  •  C  —  Um  ‘  dA 
*r—  d<y 

Dividing  the  sum  of  the  fluxes  through  the  six  surfaces 

dA\,  dAi,  dA*,  *a\,  dA'u  dA\  with  the  volume  ay  of  the  hexahedron 

we  get 


v  •  C  -  [d( 4-  a(H,H,c,n)/3x1  +  d(HiHicm)/dxt]/HlHiHji.  (1.6) 

II  —  an  (at,  Bj)cos(at,  o^)  —  nn(at,  o1)co*(nJ,  o2)  -»  un(a„  B|)cu(a2,  o2)  (1.7) 

The  projection  in  an  arbitrary  direction  m  of  the  curl  v  X  C 
of  a  vector  C  may  be  given  in  the  following  definition: 

(V  X  C)m  -  lim  M'dr,' 

*4-0  dA 

where  dA  is  the  area  enclosed  by  the  closed  boundary  of  the 
integral  h  .  Hence  it  can  be  seen  that  (v  X  C)m  is  independent 
of  the  coordinate  system.  For  example,  the  projection  of  the  curl 
of  vector  C  along  the  normal  direction  u1  of  the  surface  dA1  is 


(V  x  C)„.  -  {d[H>C, cosCnSnJJ/a*2  -  d[//2C2eo*  («*,«,)]/ 

dx*)/HiHisia  (uj,  »i) 

Finally,  we  get  the  expressions  of  the  3  components  of  the 
curl  as 


I1  ”  {3[H2C2cos(a2,  Ui)]/dx1—d[HiC)cos(ui,  a,)]/ 
dx*}/ HJHs sin  ( i*2 ,  u,)cot(ul.  Hi); 

-  {^(H^CjCos(«l,,  a2)]/a*1  -  atHtCicosCa',  at)]/ 

ajrJ}/Hs//isin(aj,  oOcosCn2,  tt2); 

I3  -  {attfjCjcosCu1,  «t)]/a*2  -  a[//2c2cos(a-\  o2)]/ 
d.r'}/,Wttf2sin  (**t,  U2)cos(o3,  a,). 


V  X  C  -  |‘n,  +  £*n2  +  §■ >us. 


(1.9) 


It  is  not  difficult  to  find  the  aero-thermodynamic  equations 
of  the  relatively  stable  flow  in  a  turbomachine  expressed  in  non- 
orthogonal  curvilinear  coordinates  by  using  vector  analysis. 

Prom  equation  (1.6)  and  v*(/>WO  —  o  ,  we  can  obtain  the  continuity 
equation  as 


d(H2H,pW'n)  djHflxpWm')  d(HxHiPW'Tl)  _ 

dxx  8xl  dxl 


(1.10) 


Substituting  equations  (1.8),  (1.9)  in  2  forms  of  the  equation 
of  motion,  i.e. 

V(  |  W\1/ 2)  -  W  X  (V  X  HO  +  2«  X  H^  —  —  -Vp/p; 

Vx(vxV)- y*(v/  -  Tvs), 

it  is  easy  to  derive  the  3  components  of  the  equations  of  motion 
12? 

u  ,  u  ,  uJ  expressed  in  terms  of  the  pressure  gradient  and  the 
gradient  of  stagnation  enthalpy  I  and  entropy  s  through  vector 
analysis: 

_1_  d(wy  +  W*  f  d[H\W\cos(u\  ai)1  _  d[HjV ,co»(u> ,  n2)l  1 
2  Hidx~  HxHi  l  3xi  dxl  ) 

_  W1  {  d[H2W2 cosCa2,  a2)1  _  d[H\W icosCa1 ,  qt)]l 
HxHi  l  dx'  dx*  ) 


+  2ulW2na(a\,  a2)cos(a3,  *,)  —  Tf^sinCu,,  ai)cos(a J,  •,)] 


/112 


(1.11a) 
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1  d(WY  j.  Wl_  f  dlH,WiCos(tt,  a2)1  _  d[HxW xcos(ux ,  gv))  j. 

T  Hidx>  HxHi  l  3*'  a*’  ] 

jy*  f  dlHilV icosCu* ,  Mi)]  _  d[HiWzcos(if,  n2)11 

h2h}1  a*2  dx3  i 

+  2u[W3sm(ui,  Ui)cos(ul,  «,)  -  ^'sinCo,,  u2)cos(a2,  «,)] 


1  d  p 

-ttVccCaz,  ^ 


1  aC#)^  J£l  I  d[HjW >cos (g* >  ifr)I  __  d\H%W xcMjtY ,  a^j 

+  a*2 

f  afW.W'.cosCa2,  aQ3  _  ajH^cosCa2,  ttiijl 
—  /fjf/,  1  dx3  dx * 


x _ i  ap 

-  u>lr  cos(B„  «r)  —  pWj  ax,  • 


(1.11b) 


(l.llc) 


T1„  |  aff/iFicosCa1,  a!)]  _  d[HiV'CQt(i?,  a,)] 

+  ^2  f  dlHjVicosCtf,  a2)]  _  d[HiVicos(u\  ai)]  | 
l  3xl  dx2  V 

Max*  axW’ 


(1.12a) 


rri  f  atHi^jCosCa2,  a2)1  _  9[Hi7|Cos(at,  at>]  | 

I  ax1  a*2  )/ 

4-  w 3  [  at/Z^iCOiCa3,  a,)]  __  a[ff2F2 cosQ2,  a2)]  | 

(  a*2  3x*  J/ 


(1.12b) 


__  wi  f  afH.FtCwCa2,  aQ]  _  d[HtV1  cotCa2,  a2)]  j ^Wj 


+  H'1' 


i^icoiCa1,  a, 

a*» 


.  /a/  79/ \ 

mmJg\dx3  dx3)' 


(1.12c) 
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2.  POTENTIAL  FUNCTION  EQUATION  OF  TURBOMACHINE  3  DIMENSIONAL 
POTENTIAL  FLOW 


Let  the  flow  in  the  interior  of  the  turbomachine  be  isentropic 
with  both  the  total  enthalpy  and  entropy  uniform  at  the  inlet. 
Hence,  (except  for  Vffv  x  V  )  vxV—  o,  $,  there  exists  a 
potential  function  ,  in  the  absolute  flow  field.  Its  gradient 
is  related  to  the  non-orthogonal  velocity  component  by 


■Hidx'  )/ 


cos  (a'  ,  if,)  (<  =  1,2,  3); 


(2.1) 


V  = 


cos(b*,  tr)  dtp 


1 _ f  cos (tt1 ,  a})  d<P 

if* ,  a,*)  I  cos  (ul,  Bj)  H^dx1  cos  (tr,  u>)  H2dx‘ 
cos(b*,  a*)  dtp  j 
cos(b3,  «i)  f/jdjru 


cos(u' 
+ 


O'  -  1.2,  3). 


(2.2) 


Substituting  the  above  equation  into  equation  (1.10),  we  get, 
after  simplification,  the  potential  function  equation 

Bu  (d'tp/dx*)  +  Bu(dl<P/dxu)  +  Bui&tp/dx*)  +  IB^d'tp/dx'dx1) 

+  2B13(Qt<P/dxldxi )  +  2Bli(dt<P/dxidxl)  +  (dtp/dxx)(dBn/dx' 

+  dBjdx1  +  3Bn/dx >)  +  (d<P/dxl)(dBu/dx'  +  dBJdx1  +  dBjdx') 

+  ( dtP/dx'XdBJdx ‘  +  dBtJdx*  +  dBJdx')  =  u>(d.h/dx' 

+  a  At/ Ox2  +  dAjdx1)  -  HMHJlOV'dlnp/Hxdx' 

+  Wld In p/ Hidx1  -1-  llr3dlap/Hidxi). 


(2.3) 


The  coefficients  Ax,  Alt  A ,,  flu.  Bu .  are  functions  of 

the  geometrical  parameters.  Its  expression  is 

„  -  cosCb'.b.),  ,  .  ncos(a‘.af'i 

Hi  co s(a?,at)  HfHj  cos(b4,  b<)cos(b',  b<)" 


(2.4) 


The  equation  for  the  potential  function  belongs  to  the  solu¬ 
tion  of  the  boundary  problems  of  the  second  kind.  Hence,  its 
boundary  conditions  may  be  obtained  from  the  condition  that  the 
product  of  the  normal  unit  vector  to  the  wall  and  the  velocity 
vector  vanishes,  i.e.,  n»W—  o  .  For  regions  up  or  down 
stream  far  away  from  the  blade,  the  direction  of  air  flow  velocity 


is  usually  parallel  to  the  inner  and  outer  shells,  e.g.  for  an 
axial  flow  turbomachine  outlet,  we  have 

Wr  —  d<P/dr  -  0  O  -  ±oc). 


For  some  special  cases,  the  potential  function  equation  may 

be  simplified  somewhat.  For  example:  by  choosing  the  angular  $ 

2 

coordinate  as  x  coordinate,  and  any  non-orthogonal  curvilinear 
coordinates  and  o,  in  the  meridional  plane  as  X1  and  X^,  we  have 


d2<P  |  H,Hi  d2<P  |  HxH,  d'Q  2  H,  cosfl»  d2<P 
Hi  sin  flu  dxu  Hi  sin0»  dx32  Hi  sin  dn  dxu  sin  du  dx'dx1 


+  [a(«^)/Hi  +  [a(S.)/a, 

_a(H,cosdu\ /Q  ,1  Kr_HiH>*. ng») 

\  dndu  )/  \  dx'  Hidx' 


-  HiHjHinn  0U 


Wx  dlnp 
Hx  dx 1 


FT,  9lnp  _F^  dlnp\ 
Hx  dx 1  Hx  dx 5  /* 


(2.5) 


The  boundary  conditions  near  the  wall  is 

n  •  W  -  n\dQ/Hxdx l)  4-  n^dt/H^x3  -  «r)  +  n\d<P/Hxdx')  -  0  (2.6) 

where  are,  respectively,  the  components  along  the 

3  directions  u1,  u2,  u^  of  the  normal  unit  vector  to  the  wall. 

3.  NON-ORTHOGONAL  CURVILINEAR  COORDINATE  STREAM  FUNCTION  EQUATION 
OF  BINARY  RELATIVE  FLOW  SURFACE 

As  early  as  1952,  Professor  Wm.  Zhonghua  proposed  a  general 
computational  model  to  obtain  the  solution  to  3  dimensional  flow 
by  taking  advantage  of  binary  relative  flow  surfaces  (S1,  S2)  to 
simplify  3  dimensional  flow  into  two  2  dimensional  problems,  so 
that  the  binary  relative  flow  surfaces  may  be  iterated.  The  flow 
matrix  method  or  velocity  gradient  method  (streamline  curvature 
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method)  may  be  used  to  solve  for  this  binary  relative  flow 
surface.  In  the  former  method,  the  stream  function  ¥  is  used 
as  the  unknown  to  find  the  various  parameters  in  the  flow  field. 

In  the  calculation  of  the  S-,  stream  surface,  it  is  common  to 

x  p 

apply  the  equation  of  motion  in  the  second  form.  Here,  x  and 
xJ  are  the  curvilinear  coordinates  in  the  stream  surface  S1, 
and  x1  is  normal  to  the  stream  surface;  the  velocity  component 
along  the  u^^  direction  and  the  curl  of  V  are  zero. 

j  u,)  ]  _  «,)![  „,) 

l  )/  (3.i) 

+  2(0  cos  (a1,  e,)/ cos  (a1,  Ui)  =  0. 

In  the  case  of  uniformly  distributed  stagnant  rotation  enthalpy 
and  entropy,  for  the  stream  surface  S1  in  the  form  of  the  surface 
of  revolution,  the  stream  function  V  that  satisfies  the  continuity 
equation 


d(HitpWl  sin  d2 
dx1 


+  'X.HjepW*  wxdg') 
dx 5 


(3.2) 


has  the  following  relations 


“  H,t pgW 2  sin  02, ,  —  —H^pglV1  sin  9 a. 

dx 1  Oxs 


(3.3) 


Substituting  the  above  equations  into  eg  (3.1),  we  may  obtain  the 
same  stream  function  equation  (Wu's  equation)  as  in  Reference  [2]: 


1  d*V  D  d2V  t,  1  d'V 

Hi  aCx2)1  h,h}  axJax5  Hi  d^y 


eJ5Li+fJ>HL„c. 

H2dxl  Hidx3 


(3.4) 


— 2  cos  0a » 

a  f 

1  H’  A 

J  ainr 

Hi  sin  0a' 

Hidx2 

a  ( 

la  Hl  \ 

ainf 

Hj0*5 

v  Wj  sin  0a' 

1  H,axJ 

d\nf  _j_ 
H,ax5 

ainf 


G  “  2uxepg  sin 


r  Hi  sin  Pa/  H>0x  **“’'"  * 

lurtpg  sin<7sin202)  +  ~2  ~  cos  9°  7*^0 

•******/{■££)• 


(3.5) 


i 
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!  The  coordinate  system  for  the  S2  stream  surface  stream  function 

1  1 

equation  may  be  chosen  as  follows:  — <p  ,  X  and  XJ  are  the  non- 
orthogonal  curcilinear  coordinates  on  the  meridional  plane.  Also 
let  n  be  the  unit  normal  vector  perpendicular  to  the  S2  stream  sur- 

^  face.  The  partial  derivatives  along  the  stream  surface  with  respect 

to  X1  and  X^  are 

„  d(  )  +  «icos(u‘,  u,)  d(  )  .  d(  )  _  d(  )  ^  n,cos(u\u3)  d(  ) 

Hidxx  Hidxx  rd(p  ’  H}dxi  Hidxi  nt  rd<p  '  (3*6) 

I 

«*  Taking  into  consideration  the  force  exerted  on  the  stream  surface 

due  to  circumferential  pressure  gradient 

/•—  —  TJL  i_  „  «  plU i  +  FiUi  +  p}U> 

J  -  and  the  perpendicular  relationship  between  the  force  P  and  the 

stream  velocity  W  as  well  as  the  relationship  between  the  unit 
normal  vector  n  and  the  stream  velocity  W,  after  simplification 
we  get  the  three  component  equations  of  the  equation  of  motion  for 
the  S2  stream  surface  as 


wx  &EL  +  w1 


^  +  >.ML\cosdn 

H,9xx  \  H\9xx  H,9xV 


H\9xx 

—  ( 2talV9  4-  «2r  4-  co* 

1  9b1 

+  w\w'  +  ^ 

-  W +  +  ft-. 

tit i  9W9  (W*.  4.  2MWlcos0t  4-  W^sinflj)  *■  F* 

H\Oxx  H&x'  \  r  ) 


(3.7a) 

(3.7b) 


or 


Wx  9(V*r)  W'  d(V9r ) 
r  Hx9xx  r  H0x' 

Wx  -lEL  +  w  +  /  ,  9WI  +  Wi  9W l) 

Hx9xx  H,9x3  V  H\9x'  H39xx) 

—  (lu>W9  4-  t olr  4-  2^sin0,  4-  tV'CW'cosdn  4-  Wl)  — 

'  r  /  H ,  Hv9x' 


-W'(iv'  +  W'cose»)±-l£L  -  Ov'ySmd3l-2%L 

H\  Ht.ax*  H39x' 


—  -  1  9f 


P  Hs9x 


7  +  f  ?• 


(3.7b' 


ny 


(3.7c) 


By  using  the  same  method,  we  may  obtain  the  equation  of  motion 
expressed  in  terms  of  the  gradients  of  the  stagnant  rotation  en¬ 
thalpy  and  entropy. 
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(3.8b) 
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\  H,9 X5  W,0*‘  )  5  v  W,dr>  / 
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(3.8c)  /116 


F*  ”  — FisinSu  - - —  — sin0j»; 

F  roqs 

c*  M  _  V  _  _  §£  a  p  • 

—  r,, 

F  rO(p 

1  9/>  «, 

F,*  —  —  F,sin0,i  —  —  F  r9(p  Bi  «n®n 
F*  -  F?  —  sin0sr,  F,*  -  F*  *■  sia0„. 

flj  "2 

Ff  W"  +  F*W9  +  Ff\Vl  -  0. 


(3.9) 


(3.10) 

(3.11) 


If  the  partial  derivatives  of  the  stream  function  of  the 
continuity  equation  which  satisfies  the  S2  stream  surface 


9{HjtpW'  sin  flu 

9xl 


g(//,rptr»Mn0„) 

9x}  0 


(3.12) 
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-  HifpfW,,sin  03i ; 
ax' 


(7* 


— Wj*7>gW,,‘sin  ®Jt. 


(3.13) 


a*3 


Substitution  of  the  above  equation  into  the  first  component 
equation  (3. 7a)  of  the  equation  of  motion  gives  the  stream  function 
equation  of  the  S2  stream  surface  [2] 


i  dlv  p  aw  i  aw  £  dv  t  F  av  _  c 
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where  the  coefficients  are  respectively 
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(3.15) 


See  Reference  [2]  for  solutions  to  this  problem  and  the  reverse 
problem. 


4.  THE  FUNDAMENTAL  SET  OF  EQUATIONS  OF  THE  Ss  AND  So  STREAM 
SURFACES  IS  SUTIBALE  FOR  THE  VELOCITY  GRADIENT  METHOD 

Professor  Wu  Zhonghua  proposed  as  early  as  1950  the  velocity 
gradient  method  (or  also  known  as  the  stream  line  curvature  method) 
[7]  with  station  by  station  calculation  to  simplify  the  2  dimension' 
al  problem  of  the  whole  field  into  the  1  dimensional  method  of 
each  calculation  station.  The  equation  for  stream  surface  S^ 
derived  in  this  paper  may  be  applied  to  the  case  of  a  surface  of 
revolution. 


The  velocity  gradient  equation  of  the  stream  surface  is 

usually  written  in  3  forms.  In  order  to  save  space,  we  shall 

write  the  velocity  gradient  equation  in  the  following  form: 

dw'/rdcp  +  PWl  4-  p  —  0  •  (4.1) 

Table  1.  Computational  formulae  for  the  coefficients  P,  Q 
of  the  velocity  gradient  equation  of  Si  and  S2  stream  surfaces 


rahsr,  ["  ttsf  <■>  *■■*«> + (~  •>  ttj 


9A  \ 


r  wr  3(Vrr) 


»A  !♦  •»’*('*' 

+  TTJ.  1  (.A  +  2CO10|,)  —  — -j-y  +  (4  +  CO»0„)  ■* 

H>*'  +  24CO10,,  +  l) 

/  9  lt>  p  ,  ,  9  In  p\\  / /  •tArm.lt..  A.  1 


fwey:  i-- coefficient ;  2~form  of  equation;  3— Si  stream  surface; 
•  4 — 32  stream  surface;  5 — Si  stream  surface. 
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V 


I 


D 


where  the  computational  formulae  for  P  and  Q  may  be  found  In 
Table  1.  In  the  table  h  —  w'/w'  represents  the  ratio  of  the 
two  velocity  components;  m'-w’/u  represents  the  Mach  number 
of  the  velocity  component  W^. 

The  Integral  form  of  the  continuity  equation.  The  flow 
quantity  of  the^stream  slab  of  thickness  t  from  9*  to  9 
is 


G  -  [9  pgxW’undnrdip.  (4.2) 

The  velocity  gradient  equation  of  the  stream  surface  S2  is 

9W'/Hx0xl  +  pw*  +  0  —  o  (4.3) 

where  P  and  Q  may  be  found  in  Table  1;  A  —  w'/w'. 


The  Integral  form  of  the  continuity  equation.  The  total  flow 
quantity  of  the  stream  slab  of  thick  ness  t  (from  h  to  t)  through 
♦  direction  1=  <4.4> 
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The  fundamental  set  of  equations  of  vicous  gas  flowing  in  turbo¬ 
machines  expressed  in  non-orthogonal  curvilinear  coordinates  will 
be  published  in  another  paper. 
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A  GENERAL  THEORY  OP  THREE-DIMENSIONAL  FLOW  IN 
TRANSONIC  TURBOMACHINES  WITH  SHOCKS* 

Xu  Jian-zhong 
(Academia  Sinica) 


Abstract 

In  order  to  consider  the  spatial  property  of  u  shock  in  a  turbomaehine  correctly 
and  in  accordance  with  a  three-dimensional  steady  transonio  flow,  the  basic  equations  - 
on  the  i»  and  i,  stream  surfaces  in  u  four -dimensional  “space”  ere  derived  through 
expanding  the  concept  of  the  stream  surface  to  the  unsteady  case.  These  equations  will  . 
be  valid  for  a  t  liree-dimensional  unsteady  flow,  too.  Applying  the  charaicteriatiea 
theory  to  the  equations,  the  characteristic  compatibility  relations  are  obtained,  and 
the  boundary  conditions  are  established.  These  equations  end  boundary  conditions 
with  the  proper  initial  conditions  entirely  define  the  problem  of  the  transonic  flows 
with  shocks  on  the  two  kinds  of  the  stream  surfaces.  Baaed  on  the  steps  of  a  complete 
method  of  solution  for  the  three-dimensional  flow  suggested  in  this  pai>er  and  using  s 
suitable  difference  scheme  the  flow  may  be  solved  numerically. 

1.  INTRODUCTION 

For  solving  the  subsonic  and  transonic  steady  3-dimensional 
flow  in  turbomachines,  reference  Cl]  proposed  the  S1  and  S2  binary 
relative  stream  surface  theory,  with  the  basic  equations  and 
boundary  conditions  of  the  binary  stream  surface  established  and 
the  method  of  solution  put  forth.  In  the  next  two  decades,  this 
method  was  widely  applied  and  verified.  In  recent  years  following 
the  development  of  transonic  turbomachines,  there  is  an  urgent 
need  for  calculating  the  3  dimensional  flow  with  spatial  shock 
discontinuity.  Based  on  the  integral  form  of  the  binary  stream 
surface  fundamental  equations,  it  was  pointed  out  in  Reference 
[2],  after  having  derived  the  relationship  among  the  air  stream 
parameters  in  front  of  and  behind  a  spatial  shock,  that  to  properly 
solve  this  type  of  3  dimensional  flow,  it  is  necessary  to  calculate 

*This  paper  has  been  presented  during  the  Second  National  Engineering 
Thermal  Physics  Conference  in  Hangzhon,  November,  1978. 


separately  on  the  binary  stream  surface  the  consistent  2  dimensional 
flow  and  Iterate  the  binary  stream  surface  calculation. 

To  establish  this  type  of  binary  stream  surface  steady 
solution  problem  on  the  stream  surface,  one  cannot  proceed  on  a 
specified  stream  surface  due  to  the  3-dimensionality  of  the 
spatial  shock.  It  is  necessary  to  extend  the  model  of  the  steady 
stream  surface  as  proposed  in  Reference  [1]  to  non-steady  flow 
and  to  merge  the  3  dimensional  flow*  to  the  binary  stream  surface 
in  the  4-dimensional  "space"  including  time  t.  In  this  paper,  we 
first  briefly  introduce  the  non-steady  stream  surface  in  this  4- 
dimensional  "space"  (or  non-constant  stream  surface)  and  point  out 
its  relationship  with  each  of  the  corresponding  instantaneous  3- 
dimensional  stream  surfaces. 

Then,  in  the  2nd  and  3rd  parts  of  this  paper,  the  basic  equations 
on  the  binary  stream  surface  consistent  with  the  3-dimensional 
flow  are  established.  They  are  suitable  for  both  non-steady  flow 
and  steady  flow  with  spatial  shock.  In  the  4th  part  of  this 
paper,  we  derive  the  characteristic  consistency  relationship  on 
the  binary  stream  surface  based  on  characteristic  curve  theory  and 
propose  the  set  of  boundary  conditions  on  the  up  and  down  stream 
boundaries  and  on  the  object  surface  suitable  for  solving  steady 
state  motion.  The  initial  conditions  may  be  given  according  to  the 
usual  method.  Thus,  the  steady  solution  problem  has  been  completely 
established. 

Lastly,  we  briefly  describe  the  procedures  to  solve  3-dimensional 
flow  by  iterating  the  binary  non-steady  stream  surface  calculations. 

We  would  like  to  express  our  gratitude  to  Professor  Wu  Zhonghua  for 
his  guidance  in  developing  this  paper. 


It  has  been  pointed  out  by  Huang  Ruixin  that  the  steady  stream 
surface  may  be  extended  to  non-steady  stream  surface. 


2.  THE  NON-STEADY  STREAM  SURFACE  AND  ITS  GEOMETRICAL  /121 

RELATIONS 


A  spatial  stream  surface  may  consist  of  the  collection  of  all 
stream  lines  through  the  points  on  a  curve  that  is  not  a  stream 
line  in  a  flow  field. 


In  non-steady  flows,  the  stream  surface  changes  with  time. 
In  this  paper,  we  shall  describe  this  non-steady  stream  surface 
in  the  4  dimensional  "space”  including  time  t:  }(r,  <p,  *,  <)  —  0. 


Thence , 


dr  dqf»  da  dt 


0 


and  its  unit  normal  vector  N(Nr%Nm%N..>N,)  satisfies 


Thus,  on  the  stream  surface  in  the  4-dimensional  "space", 

NJr  +  N.rdtp  +  N,Ja  +  N,Uit  —  0, 


N, _ 1  dr  N, _ r  (!) 

N,  U  Bt*  N.  U  Bt 

where  the  dash  line  in  a  derivative  indicates  the  partial 
derivative  along  the  stream  surface. 

On  the  other  hand,  there  is  a  stream  surface  *(r* *»*) “ °» 
at  every  instant  in  3-dimensional  space  with  its  unit  normal  vector 

••)-  as  a  function  of  time  t.  At  the  same  time,  according 

to  the  definition  of  the  stream  surface,  we  have  n'W  —  d,  at 
every  instant. 

It  should  be  pointed  out  that  the  following  relations  hold 
between  the  unit  normal  vector  of  the  3-dimensional  spatial  stream 
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surface  and  the  4-dimensional  " spatial”  stream  surface  n  and  N: 


NJN,  •  njn ,,  NJN ,  — •  njn,,  NJN.  -■  nmJnm  (  2  ) 

With  the  above  relations  for  a  non-steady  stream  surface,  we 
can  then  merge  the  basic  set  of  equations  of  absolute  thermal 
motion  for  a  non-viscous  gas  expressed  in  a  relative  coordinate 
system  rotating  with  constant  angular  velocity  u  given  by 
Reference  [1] 
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dt  r  dr  r  dtp  d t 


aw,  +  W'  aw,  +  w.  awj  +  w 

dt  'dr  r  dtp 


0 


_  J_  Pjl 

P  dr 


aw, 

t  +  SEi 

>  4. 

Wm 

aw. 

,  4- 

W 

aw. 

wrw 

i  4.  2 taW 

at 

'^T 

■  T 

r 

dtp 

*  i 

**  a 

a. 

r 

_ _ 1_  dp 

pr  dip 

aw 

I-  +  W  Mi 

L  4- 

w . 

aw. 

,  4. 

W 

aw.  — 

_ I_ 

$L 

dt 

» 

,  r 

dtp 

*  i 

**  a 

dt 

P 

dt 

dl_ 

wm 

dl_ 

+  W. 

at_ 

0 

dt 

dr 

r 

dtp 

a * 

p  at 

(3) - 

(4) 


(5) 

(6) 

(7) 


onto  the  non-steady  stream  surface  in  binary  4-dimensional  "space" . 
It  should  be  explained  that  as  a  closed  system  by  equations,  besides 
the  above  equations,  the  equation  of  state  for  a  gas 

p-pRT  (8) 


and  the  relations 


dh~CJT,  imk  +  W'/l-U'/l  (9) 


should  be  Included. 

Since  these  algebraic  equations  retain  their  forms  on  the 
stream  surfaces,  we  will  not  specifically  describe  them  in  the 
following  sections. 


3.  BASIC  EQUATIONS  ON  THE  NON-STEADY  STREAM  SURFACE  S1 

By  extending  the  definition  of  the  stream  surface  St  in  [1] 
to  the  4-dimensional  "space”,  it  may  be  expressed  as  ^  (r,<p  *  t)  ■ 
Any  function  q  on  this  surface  may  be  expressed  as 

q  —  *»  O  —  t,  l),<p,  t,  t  ]. 


We  notice  that  in  relation  (2),  the  partial  derivative  of  q  on 
the  stream  surface  S.^  may  be  written  as 


1  dq  _  1  dq  _ dq  3q  _  dq  _  ffj  dq  3q  _  _  l/A/t 

r  dtp  r  dtp  n,  dr  '  3s  ds  n,  dr  *  3i  dt  Nr  dr 


(10) 


Int roducing 


2ln»  \\n.diW.r)  dW  .  .  dWA  ^  UfjU  /Mag.  +*±*1 

~T  TrV7'^r^  +  ■"  Or  +  dr  ]  V,  dr  r  (11) 

the  continuity  equation  (3)  may  be  written  as 

SQ>p)  t  1  5(bpW.)  +l_  BjrbpW^  _  0  (12) 

dt  r  dq>  r  ds 

Similarly,  equations  ( 4 ) — ( 7 )  may  be  expressed  respectively  as 


3 w.  (  wm  dw, 

3t  r  dtp 


f.r  B\Vr  K* 

Wmir*  ?  * 


3tvm  t  y,  3wm  |  lt,  fly.  tVjy. 
3<  r  *  3z  r 


+  2wWr  — 


fly.  ^  y. 

9(  r  dtp 
3i  ^  y.  di 
3t  r  dip 


pr  dtp 


(13) 


(14) 


(15) 


(16) 


where 


r--  -L*Ln-H*LW  ,  _  uat,  /  a/  j_ 

*rP  dr  iV,  dr  *  ivr  Vdr  P  dr ) 


(17) 


If  we  consider  steady  flow,  5/fli  —  o, N,  —  o,  ,  then  the 
above  equations  become  respectively  the  basic  equations 
(34b)*  (40),  (41),  (39*)>  (39b)  and  (39c).  on  the  stream  surface 


in  Reference  [1].  In  other  words,  compared  to  the  equations 
on  a  steady  stream  surface  S^,  besides  the  extra  term 

in  the  non-steady  stream  surface  equations,  time  related  terms 

UN,  din  ft  UNi 

and  ~N/~ij^~*  are  included  in  the  expression 
of  2\*b/dt  and  ^  and  terms  like  J  appear  in  the  energy  equation. 


We  see  that,  different  from  solving  the  basic  equations  of 
transonic  flow  with  shock  discontinuity  on  S1  stream  surface, 
in  the  equations  (12)-(l6),  b  includes  non-steady  terms  which 
change  with  time  in  the  form  of  ^  ;  at  the  same  time  in  f  and 
J  there  are  effects  of  non-steady  factors.  These  are  results  due 
to  the  adoption  of  the  model  of  a  non-steady  stream  surface  in  4- 
dimensional  "space”  when  we  derive  these  equations  by  treating 
the  3-dimensionality  of  shock  discontinuity  in  this  paper. 


Taking  the  fact  into  consideration  that  in  the  process  of 
numerical  solutions,  it  is  more  accurate  to  use  the  difference 
scheme  in  conservative  forms  through  the  shock  discontinuity,  we 
transform  equations  (12)-(l6)  into  the  following  weak  divergence 
form  (corresponding  equations  may  be  obtained  by  using  I  and  S) 

9i  dtp  5*  iV,  ( 1  o 


as=i+ +ir('*  7) +*l-*  ■  “ '(7 'A 


(19) 


lurm 


*  -  n-p\V.,l- pWt. 


(21) 

(22) 

(23 )— (25 ) 
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4.  BASIC  EQUATIONS  ON  THE  NON-STEADY  STREAM  SURFACE  S2 


By  extending  the  steady  ,S2  stream  surface  in  3-dimensional 
space  to  the  4-dimensional  "space”,  it  may  be  written  as 

h(r>  <p>  *»  *)  —  b  •  On  this  surface  any  function  q  may  be 

expressed  as  <?  —  *('.  <p»  *»  0  tp(r,  *,  t),  *,  <],  and 

the  partical  derivative  on  the  stream  surface  may  be  transformed 
into 


On.  «.  _  jiL  d<?  _  »*  d« 

S'”  dr  a'.r  dip  ’  3*  dc  ftp  * 


Otf  —  dg  _  owV'  9^ 

3<  d/  AT.  ~dj> 


(26) 


Introducing 


i//  n.r  \  dtp  dip  dtp  /  N"m  dtp  r 

we  may  write  the  continuity  equation  (3)  as 

+  K>'pu\)  dp'riv.) 

dt  dr  3, 

Similarly,  equations  ( 4 ) — ( 7 )  are  respectively 


(27) 


Hr- 
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dr  9,  |9  3x 

5/  *3*  ^'3r+V’37  r 


(28) 


(29) 

(30) 
(3D 
(32) 


where 


F-  -J-±  d£.n'  _  ja^qjV  ■  r__  wtT.  (  9»  >  d,  \ 

*»r  P  dip  tfu  dip  ’  AT.  'd  ~P  dip  ' 


(33) 


Just  as  the  case  of  the  S^  stream  surface  in  steady  flow 
the  above  equations  are  transformed  into  equations  (100),  (98), 


/12 
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(99a),  (96a),  (96b)  and  (96c)  in  Reference  [1];  and  the  difference 
between  the  steady  S2  stream  surface  and  the  non-steady  Sg  stream 
surface  is  also  similar  '.0  the  case  of  the  stream  surface. 

The  continuity  equation  (28)  is  already  at  the  strong  divergence 
type.  The  other  equations  may  be  transformed  into  the  following 
weak  divergence  form  (also  expressible  in  terms  of  I  and  S): 
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We  should  point  out  that  the  basic  equations  on  a  binary  non¬ 
steady  stream  surface  above  are  suitable  for  both  solving  the  non¬ 
steady  flow  and  the  steady  flow  with  shock.  The  difference 
in  the  solutions  of  these  2  problems  is  in  different  boundary 
conditions . 


5.  BOUNDARY  CONDITIONS  AND  INITIAL  CONDITIONS 

In  the  previous  2  sections,  we  have  established  the  hyper¬ 
bolic  partial  differential  equations  on  the  binary  non-steady 
stream  surface.  In  the  following,  the  corresponding  boundary  and 
initial  conditions  for  solving  steady  motion  are  given.  In 
initial  value  —  boundary  value  computations,  the  parameters, 
on  the  object  surface  and  the  up  and  down  stream  boundaries 
are  related  to  those  in  the  flow  field  at  a  previous  instant. 
Thus,  one  should  derive  the  characteristic  consistency  relations 
according  to  the  characteristic  curve  theory  [3]  of  hyperbolic 
equations  and  then  determine  the  appropriate  boundary  conditions 
[4-71. 
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(1)  Boundary  conditions  on  the  S.^  stream  surface. 


The  following  consistency  relations  may  be  obtained  on  the 
S]  stream  surface: 
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where  y  represents  the  angle  between  the  z  axis  and  the  exterior 
normal  at  a  point  on  the  boundary.  Equations  (38)  and  (39)  are 
the  flow  characteristic  consistency  relations,  while  (40)  and 
(4l)  are  respectively  the  wave  characteristic  consistency  relations 
of  propagation  with  velocity  -a  and  "a”  relative  to  the  gas. 


For  the  upstream  boundary:  When  the  meridional  plane 
velocity  W,A<»  ,  only  (40)  may  be  used.  Thus  3  parameters  need 

be  specified  and  the  other  parameter  should  be  derived  from  (40). 
Notice  that  now  /*“•*  ,  and  the  3  specified  parameters  (e.g.  , 

T,  and  air  stream  angle)  do  not  vary  with  time  while  the 
parameters  are  all  uniform  along  the  9  direction.  Equation  (40) 
will  be  greatly  simplified. 

For  the  down  stream  boundary:  Here,  if  the  outlet  air  stream 
is  subsonic,  then  eg.  (38),  (39)  and  (41)  are  usable  consistency 
relations.  We  need  to  specify  another  parameter  (i.e.  the  reverse 
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pressure  p).  Notice  that  — 0  and  the  3  equations  above  may 
be  greatly  simplified.  If  the  outlet  air  stream  is  transonic, 
the  4  consistency  relation  can  all  be  used  and  no  parameters  can 
be  specified. 

For  the  blade  surface:  The  relations  (38),  (39)  and  (41) 
can  all  be  used.  The  specified  condition  is  that  the  stream 
velocity  should  be  tangential  to  the  surface. 

Besides,  there  is  still  the  periodic  condition  for  the 
Si  stream  surface  calculation  which  may  proceed  in  a  way  similar 
to  an  interior  point.  Of  course,  when  the  shock  extends  beyond 
the  passage  way,  we  need  to  analyse  more  concretely  whether  to 
choose  periodic  conditions  for  pressure  surface  or  for  suction 
surface. 


(2)  Boundary  conditions  on  the  S2  stream  surface. 


Similarly  the  consistency  relations  now  are 
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where  represents  the  angle  between  the  z  axis  and  the  exterior 

normal  at  a  point  on  the  boundary.  Equations  (42)  and  (43)  are  the 


flow  characteristic  consistency  relations  while  (44)  and  (45) 
are,  respectively,  the  wave  characteristic  consistency  relations 
with  propagation  velocity  -a  and  "a”  relative  to  the  gas. 

For  up  stream  boundary:  When  ,  only  (44)  applies. 

Thus,  3  parameters  need  be  specified  with  the  other  parameter 
derived  from  (44).  Notice  that  now  m'  and  the  3  specified 
parameters  (e.g.  h* and  M  number)  do  not  vary  with  time,  while 
the  parameters  along  the  4  direction  are  all  uniform.  E.g.  (44) 
will  be  greatly  simplified. 

For  down  stream  boundary:  Now  if  the  outlet  air  stream  is  sub¬ 
sonic,  then  (42),  (43)  and  (45)  are  usable  consistency  relations. 

One  parameter  needs  by  specified  (e.g.  inverse  pressure  p  or  vmr )  ). 
Notice  that  now  **'  “  0  and  the  above  3  equations  will  also  be 
greatly  simplified.  If  the  outlet  air  stream  is  transonic,  the 
4  consistency  relations  can  all  be  used  and  no  parameters  can  be 
specified. 

For  wall  surface:  (42),  (43)  and  (45)  all  apply.  The  /126 

specified  boundary  condition  is  that  the  air  stream  velocity 
should  be  tangential  to  the  wall  surface. 

(3)  Initial  conditions 

It  has  been  proved  In  theory  [8]  that  the  choice  of  initial 
conditions  does  not  affect  the  final  steady  state  solution.  But, 
if  the  choice  is  not  proper  at  all,  then  an  excessive  oscillation 
may  be  induced  and  the  stability  of  the  numerical  solution  will  be 
destroyed.  Also,  the  choice  of  initial  conditions  will  affect 
the  speed  of  reaching  stable  solutions.  Taking  these  factors  into 
consideration,  one  may  choose  the  subsonic  solution  or  the  para¬ 
meter  distribution  that  satisfies  the  continuity  equation  on  the 
geometrical  center  line  while  smoothly  varying  elsewhere  as  the 
initial  condition. 
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i. 


PROCEDURES  FOR  A . COMPLETE  SOLUTION  OF  3-DIMENSIONAL  FLOW 


In  the  above  we  have  transformed  the  problem  of  solving  steady 
state  3-dimensional  flow  with  shocks  in  turbomachines  into  the  initial 
boundary  value  problem  for  2  families  of  3  variable  hyperbolic 
equations,  and  suggested  the  way  to  specify  initial  and  boundary 
conditions.  From  these  we  have  completely  established  the  problem 
of  finding  solutions  on  the  binary  non-steady  stream  surface  in 
4-dimensional  "space".  To  obtain  the  solution  to  the  3-dimensional 
problem,  besides  obtaining  respectively  the  stable  solutions  of  the 
initial  value-boundary  value  problem  on  the  2  families  of  stream 
surfaces,  it  is  also  necessary  to  Iterate  them  until  convergent. 

To  save  time,  the  solution  process  may  adopt  the  method  of 
both  time  step  and  stream  surface  iteration:  After  each  time  step 
a or  several  steps  2a**  ,  an  iteration  of  the  2  stream 

surface  families  is  carried  out,  i.e.,  taking  the  result  of  each 
stream  surface  calculation  and  the  shape  of  the  4-dimensional  stream 
surface  [e.g.  (1)]  and  using  them  in  the  calculation  of  the  next 
time  step  length  for  the  other  stream  surfaces.  This  is  re-cycled 
whether  the  solution  in  each  calculation  has  reached  stability 
or  not.  When  a  stable  solution  is  approached  at  each  time  step 
length  one  must  iterate  the  stream  surfaces.  Finally,  when  a 
stable  solution  is  obtained  on  each  stream  surface,  the  final 
solution  for  the  whole  3-dimensional  flow  is  found.  This  solution 
procedure  clearly  shows  that  transonic  flow  with  shock  cannot 
be  solved  on  a  specified  stream  surface  as  pointed  out  in  Reference 
[2].  In  other  words,  it  is  precisely  through  the  continuous 
adjustment  of  the  stream  surface  position  and  shape  during  the 
solution  process  that  one  can  obtain  the  2-dimensional  shock  in 
agreement  with  the  spatial  shock  on  the  binary  stream  surface. 

The  above  solution  procedure  closely  correlates  the  calculated 
shock  on  a  single  stream  surface  and  the  spatial  approximation 
of  the  shock  which  shortens  the  time  required  to  reach  a  stable 
solution. 


7.  CONCLUSIONS 


Based  on  the  extension  of  the  stable  stream  surface  model 
in  3-dimensional  space  to  the  unstable  stream  surface  in  4- 
dimensional  "space”,  in  this  paper  we  transformed  the  partial 
differential  equations  of  non-viscous  gas  absolute  thermal  motion 
in  turbomachines  into  the  binary  non-steady  stream  surfaces  3^,  and 
S2,  and  established  the  corresponding  basic  equations.  Compared 
to  the  equations  on  stable  stream  surfaces,  in  these  equations, 
there  is  a  contribution  to  b  (or  b')  and  "stream  slab  force"  due 
to  the  shape  of  non-steady  stream  surface,  besides  the  partial 
derivatives  which  include  the  thickness  of  the  stream  slab  in 
4-dimensional  "space"..  In  the  energy  equation  related  terms  J 
(or  J^)  also  appear.  These  equations  are  suitable  for  both  3- 
dimensional  non-steady  flow  and  steady  flow  with  shocks. 

It  should  be  pointed  out  that  the  equations  for  solving 
transonic  flow  on  stream  surfaces  (mainly  surface  of  revolution) 
in  both  domestic  and  foreign  references  are  all  established  on 
steady  stream  surfaces  that  do  not  vary  with  time.  The  thick¬ 
ness  of  the  stream  slab  is  also  independent  of  time.  Hence,  the 
shock  thus  calculated  is  not  part  of  the  real  spatial  shock. 

In  fact,  the  basic  equations  in  these  references  are  not  obtained 
from  the  basic  equations  for  solving  spatial  shocks.  Their  results 
fail  to  consider  the  close  relationship  between  stream  surfaces 
of  the  same  family  and  between  the  stream  surfaces  of  the  2  families 
induced  by  the  3  dimensionality  of  the  shock. 

According  to  characteristic  theory,  we  have  derived  the 
characteristic  consistency  relation  on  binary  stream  surfaces  and 
.the  boundary  conditions  for  steady  flow.  After  suitably  choosing 
the  initial  conditions,  the  problem  of  finding  solutions  for  trans¬ 
sonic  flow  with  shocks  on  binary  stream  surfaces  is  completely 
established. 
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Taking  into  consideration  the  tight  correlation  between 
stream  surfaces  with  shocks,  for  obtaining  the  solution  of  3- 
dimensional  flow,  we  may  use  the  method  of  simultaneous  time  step 
and  stream  surface  iteration.  It  can  save  computation  time. 

Thus,  one  may  find  the  numerical  solution  after  suitably  selecting 
the  difference  scheme. 
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AN  EXPERIMENT  TO  IMPROVE  THE  SURGE  MARGIN  BY 
USE  OF  CASCADE  WITH  SPLITTER  BLADES 

Tsui  Chih-ya,  Zhou  Sheng,  Lu  Ren-fu,  Zhang  Lian-chi 
(Beijing  Institute  of  Aeronautics  and  Astronautics) 


Abstract 

a 

Through  experimental  comparison  cn  a  basis  of  equal  profit^  sectional  area  p**r 
unit  frontal  length,  which  implies  a  basis  of  equal  weight,  a  cascade  with  splitter 
blades  is  found  to  give  a  maximum  deflection  3.5  degrees  greater  than  an  ordinary 
cascade  tested.  The  corresponding  incidence  angle  is  6.5  degrees  greater,  thus  mar¬ 
kedly  improving  the  surge  margin.  The  nominal  deflection  and  thus  work  addition 
ability  also  increase. 

Prediction  is  also  made  on  a  basis  of  equal  solidity  that  the  use  of  splitter  blades 
would  give  about  the  same  nominal  work  addition  ability  as  an  ordinary  cascade,  but 
weight  is  substantially  reduced  to  70%  approximately. 


1.  FORWORD 

It  is  well  known  that  the  maximum  deflection  may  be  increased 
by  increasing  the  solidity  of  the  pressure-resisting  cascade, 
hence  improving  the  surge  margin;  however,  the  weight  also 
increases  and  the  efficiency  decreases  due  to  greater  flow 
resistance.  The  basic  mechanism  is  mainly  due  to  the  fact  that 
with  a  high  solidity  and  closer  blade  separation,  the  pressure 
difference  decreases  when  the  same  air  stream  turns  and  the 
minimum  pressure  at  the  rear  side  of  the  blades  together  with 
the  reverse  pressure  gradient  is  more  gradual,  so  that  separation 
will  only  occur  at  higher  attack  angle  and  rotation  angle,  hence 
the  increase  in  maximum  angle  of  rotation  and  surge  margin. 

It  often  happens  that  the  minimum  pressure  point  of  a  low 
speed  cascade  blade  back  shifts  forward.  This  offers  a  clue. 

But,  then  why  should  the  solidity  be  increased  along  the  entire 

9 
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blade  chord? 


In  this  investigation,  we  want  to  study:  for  a  cascade 
with  splitter  blades  which  has  the  same  profile  sectional  area 
per  unit  frontal  length  as  an  ordinary  cascade  (Figure  1)  i.e. 


Figure  1.  Ordinary  cascade 
and  cascade  with  splitting 
blades. 

Key:  1 — ordinary  cascade; 

2 — cascade  with  splitting 
blades. 


Figure  2.  The  attack  angle 
characteristics  of  2  types 
of  cascade. 

Key:  1 — A  ordinary  cascade; 

2 — o  cascade  with  splitting 
blades;  3 — +  angle  of  rotation. 


F«  +  F*  _  F* 

t'  t 
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whether  the  maximum  angle  of  rotation  and  the  corresponding  angle 
of  attack  will  increase  for  the  same  weight  of  the  blades. 
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2 .  EXPERIMENTAL  SET-UP 


The  experimental  horizontal  tyoe  cascade  test  outlet  120  X 
490,  and  front  blade  .w,-o.3,  (2.3-4)xio*. 

The  middle  section  40mm  blade  front  and  rear  total  static  pressure 
non-uniformity  is  less  than  1.5?  and  the  non-uniformity  of  the 
blade  front  and  rear  air  stream  angle  is  less  than  1°.  The 
periodicity  of  the  3  central  channels  along  the  cascade  distance 
is  good,  which  provides  reliability  for  the  comparison  study. 

The  ordinary  cascade  consists  of  10  blades  with  blade  type 
10C4/30CS0  ,  chord  length  60mm,  cascade  distance  46.2mm  (solidity 

1.30),  angle  of  Incidence  60°  and  is  made  of  glass  epoxy  with  an 
epoxy  coating. 

Two  blades  are  removed  from  the  ordinary  cascade  and  replaced 
by  7  small  131  blades  with  chord  36mm,  inserted  uniformly.  The 
distances  from  the  end  to  the  large  blades  on  the  2  sides  are  both 
31.4mm.  The  angle  between  the  blade  chord  and  the  large  blade 
chord  is  5°  which  makes  the  angle  of  attack  ismallx  liarge  "  l0* 

The  1  dimensional  divergence  angles  on  the  2-side  channels  are 
approximately  equal.  The  small  blade  is  made  of  aluminum  or 
brass . 

When  the  4-hole  sensing  head  is  adjusted  in  the  32  m/s  wind 
tunnel,  the  error  in  the  total  pressure  1.5?,  the  error  in  static 
pressure  4.3?  and  the  error  in  the  direction  0.25°.  The  positions 
of  measurement  are  respectively  at  2  chord  lengths  in  front  of 
the  cascade  and  one  chord  length  behind  the  cascade.  Measurements 
at  many  points  along  the  cascade  are  averaged.  Pressure  distri¬ 
bution  is  measured  at  the  central  cross-section  of  the  central  blade. 

3,  RESULTS  AND  DISCUSSIONS 

The  characteristics  of  the  attack  angles  of  the  2  cascades 
are  shown  in  Figure  2.  The  typical  blade  surface  pressure  distri¬ 
bution  t,  p  —  P\  is  shown  in  Figure  3  a,b,c,d. 
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Compared  to  that  of  the  ordinary  cascade,  the  maximum  angle  of 
rotation  of  the  cascade  with  splitting  blades  is  35. 3° ,  i.e.  3.5° 
more,  and  the  angle  of  attack  is  19.8°,  i.e.,  6.5°  more.  Although 
the  loss  coefficient  5  of  0.435  is  .103  larger  than  that  of  the 
ordinary  cascade  due  to  the  increase  of  frictional  surface  and 
of  wake  mixing,  yet  the  separation  is  apparently  delayed  and 
surge  margin  improved.  The  tendency  of  the  pressure  distribution 
curve  is  also  concrete  evidence.  As  pointed  out  by  the  experiment 
in  Reference  [1],  the  larger  the  axial  flow  is  than  W'uri/M'i.rj  ,  the 
greater  will  be  the  reverse  pressure  gradient  on  the  back  siie  of 

the  blade  and  the  easier  will  the  air  stream  seoarate.  Now  -rr  ■* u  ' 

\v  u 

at  maximum  angle  of  rotation  for  the  cascade  with  splitting  blades 
is  slightly  greater  than  0,745  at  maximum  angle  of  rotation 

for  the  ordinary  cascade  which  demonstrates  that  when  compared 
according  to  the  same  flow  ratio,  the  advantages  of  the  cascade 
with  splitting  blades  will  be  slightly  enhanced. 

The  nominal  angle  of  rotation  for  the  ordinary  cascade 

-  0.8  x  31.8  -  25.5°  .  Here  i*  -  2.5*.  fi?  -  42.5°,  ft*-  68°; 

From  the  curve  for  synthetic  nominal  angle  of  rotation  according 
to  this  and  solidity  1.3  (as  in  Reference  [2]),  A0*  is  also 
25.5°,  just  as  in  our  experiment. 

The  nominal  angle  of  rotation  for  the  cascade  with  splitting 
blades  49*  —  0.8  x  35.3  —  28.3°  »  which  is  2.8°  greater  than  that  of  the 

ordinary  cascade.  The  loss  coefficients,  however,  are  basically 
the  same  as  can  be  seen  from  Figure  2.  Hence,  the  work  addition 
capability  of  the  nominal  angle  of  rotation  is  also  improved. 

Another  very  meaningful  comparison  may  be  predicted.  If  we 
replace  the  small  blades  with  large  ones,  then  the  solidity  — _  —  1.9 

The  airstream  off-path  angle  may  be  found  by  using  the  empirical 
formula 


8m 


(75°  -  68°),^  -  5-*# 
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i.e.  ft  mm  75°  —  5.8°  —  69.2°  .  From  Reference  [2],  we  find 

—  29.5°  ,  close  to  that  for  the  cascade  with  splitting  blades, 

while  the  cross-section  (or  weight)  of  the  cascade  with  splitting 
blades  reduction  ratio  is  CIO*  -H-  6*):(10*  10*)  ■—  0.68. 

This  has  raised  some  new  questions:  Is  it  possible  that  by  changing 
every  other  blade  in  an  ordinary  cascade  which  has  been  in  use  for 
years  to  a  small  one,  the  work  addition  will  remain  the  same,  while 
the  weight  and  cost  will  be  much  reduced? 

The  only  experiment  independently  carried  out  in  the  sixties  is  in 
Reference  [3].  They  tried  to  use  small  blades  at  different  positions 
in  the  middle,  or  rear,  of  the  passage,  but  the  chord  length 
is  only  25$  of  that  of  the  big  blade.  The  practicality  of  technology 
and  strength  is  far  inferior  to  that  in  our  paper.  The  characteristics 
of' the  attack  angle  are  not  given,  and  the  pressure  distribution  of 
the  small  blade  is  not  measured. 

In  the  seventies.  Reference  [4,5]  experimented  with  supersonic 
extended  pressure  cascades  with  small  blades  added  from  axial  chord 
length  50$  to  rear  edge.  At  air  Intake  Mach  number  about  1.46,  the 
air  outlet  angle  is  made  to  be  8-9$  closer  to  the  axis,  but  the 
loss  is  increased  to  25$.  In  Reference  [6],  a  supersonic  compressor 
is  designed  with  small  blades  added  in  the  rear  half  of  the  passage 
to  reduce  the  air  outlet,  off-path  angle  of  the  rotor,  and  to  strive 
toward  a  single  stage,  total  pressure  ratio  of  3-0.  There  is  a  pre¬ 
liminary  report  (Reference  [7]),  but  the  intention  there,  is  apparently 
different  from  our  paper. 

All  in  all,  our  study  is  only  a  comparison  experiment  at  low 
speed  for  cascade  distance,  angle  of  bending,  angle  of  attack,  relative 
size  and  position  of  the  small  blade.  The  result  is  encouraging. 

It  is  necessary  to  further  develop  experiments  to  compare  things 
in  more  detail  so  as  to  provide  a  practical  reference  in  the  design 
for  compressors. 
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THE  PERFORMANCE  CALCULATION  OF  AN  AXIAL 
FLOW  COMPRESSOR  STAGE 


Zhang  Yu-jing 


Abstract 

Based  on  the  simplified  two-dimensional  flow  equations,  the  computer  program 
for  the  direct  problem  is  developed.  By  means  of  this  program  several  subsonic  axial 
flow  compressor  stages  are  calculated.  The  calculated  results  are  compared  with  ex¬ 
perimental  data.  This  calculation  method  is  evaluated  and  problems  to  be  solved  in 
the  performance  calculation  are  presented. 


1 .  FORWORD 

The  calculation  of  the  axial  flow  compressor  problem  means  the 
calculation  of  compressor  characteristics  when  all  the  geometrical 
parameters  of  the  compressor  are  known.  This  kind  of  calculation 
is  important  and  meaningful  whether  to  determine  the  combined 
operation  line  of  the  compressor  and  the  turbine,  or  to  understand 
the  properties  of  the  rotatable  airfoil.  If  the  characteristics 
can  be  computed  with  good  enough  accuracy,  then  less  effort  can  be 
spent  on  testing  and  adjusting  or  tuning  the  compressor.  Because 
of  our  research  on  the  attached  surface  layer,  the  surge  and  the 
properties  of  the  cascade,  etc.  are  not  enough.  Therefore,  the 
compressor  characteristics  are  still  estimated  approximately,  or 
obtained  from  experiment.  There  are  m.'.ny  problems  requiring  further 
studies  in  the  method  of  calculating  the  characteristics  of  the  whole 
machine  based  on  the  binary  relative  flow  surface  theory  [1,2]. 

There  is  still  a  considerable  distance  before  actual  engineering 
applications  of  that  method.  However,  this  is  a  general  method 
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to  understand  the  gas  flow  inside  the  compressor,  and  is  the 
direction  in  which  the  computation  of  compressor  characteristics 
moves . 

To  study  problems  in  this  area,  we  developed  a  computer 
program  in  1975  for  the  compressor  problem,  carried  out  some  cal¬ 
culations  for  some  subsonic  compressor  stages,  and  compared  the 
computational  results  to  the  current  experimental  results.  This 
was  a  preliminary  attempt  to  calculate  the  characteristics  with 
the  simplified  2  dimensional  flow  method.  Thus,  we  did  not  place 
as  much  weight  on  the  equations  themselves,  and  instead  we  con¬ 
centrated  on  such  problems  as  the  convergence  rate  in  the  calculation, 
the  effect  of  the  wrong  choice  of  some  empirical  laws  (e.g.  determina¬ 
tion  of  the  reference  attack  angle)  in  the  design  and  in  the  computa¬ 
tion,  etc.  This  made  us  realize  that,  in  order  to  obtain  reliably 
computed  characteristics,  one  must  have  complete  information  on  the 
accurate  cascade  reference  attack  angle  and  the  corresponding  angle 
of  lag,  as  well  as  the  computational  method  for  the  cascade 
characteristics  under  reasonable  work  conditions.  These  two  aspects 
are  more  important  than  an  accurate  loss  model. 

According  to  the  results  of  calculations  in  this  paper,  it  is 
feasible  to  use  the  computational  method  for  subsonic  compressor 
stage  characteristics.  Further  computations  and  experimental 
verification  are  needed  for  making  calculations  on  the  whole  compressor 


2.  BASIC  ASSUMPTIONS  AND  EQUATIONS 

In  this  calculation,  we  assume  the  flow  to  be  stable,  axial 
symmetric  and  adiabatic.  The  fluid  is  perfect;  viscosity  is  neglected 
but  the  change  in  the  entropy  due  to  it  is  taken  into  consideration. 
The  rotating  angular  speed  of  the  rotor  is  constant.  This  greatly  . 
simplifies  the  problem. 


3 


/ 

In  the  v( r,<p,2 )  coordinate  system,  the  continuity  equation 
and  equation  of  motion  are  respectively 


The  energy  equation  is: 

Flowing  through  a 

Flowing  through  a 
The  equation  of  state  is 


1  d((>wtr )  J  d(pW.)  __  0 
r  dr  dz 

(1) 

WrM£jiLL  a-  W t  — * '-C-ur'>  —  0 
dr  dz 

(2a) 

W  aw’  +  w.dEs-c±  -  -  So  dp 
'dr  ds  r  P  dr 

(2b  i 

IV  MV.  +  iv  dW-  ~-Sod£_ 
dr  dz  P  dz 

(2c) 

moving  airfoil:  Dl/Dt  ■=  0 

(3) 

stationary  airfoil:  DH^Dt  =  0 

(3'  ) 

p  —  pRT 

(4) 

According  to  the  continuity  equation,  equation  of  motion,  and  the 
equation  of  state,  we  can  derive  the  radial  equation  of  motion: 

So  $£.  _  ( 1  —  M»y  £i  _j_  Jec3  £&\  +  c\un<p  d(r  tamp) 

P  dr  \\-MlJ\r  rj  1  -Ml,  rdr 

_ <5> 

1  -  M'.t,  Dl 

This  is  e.g.  (31)  of  Reference  [33  which  has  taken  into  account 
the  entropy  change. 

To  make  the  computation  easier,  we  introduce  the  stream  line 
parameter  A,  a  »  1  +  tanJ<p  +  tan2£, 

Here  tan  fi  —  w  Jw 
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Equation  (5)  may  be  written  as  follows: 


[3/ 

A 

2W\ 

A 

l2  1 

_tan£ 

9( 

For  the  computational  cross-section  behind  a  stationary 
airfoil,  it  is  only  necessary  to  take  i  »  ”  °»  P 


Eq.  (6)  is  the  radial  equation  of  motion  which  is  to  be  solved 
in  the  characteristic  calculation. 


3.  SOLVING  THE  EQUATIONS 

We  use  the  iteration  method  to  solve  the  equations.  Some 
explanations  of  the  problem  in  the  solution  are  as  follows: 

1.  Initially  when  the  stream  line  is  assumed,  the  radial 
point  is  determined  according  to  the  equal  ring  surface. 

2.  When  '<P;  and  *m  are  calculated,  the  shape  of  the 
stream  line  is  taken  to  be  a  cubic  function.  In  each  stage  of 

the  calculation,  there  are  3  known  positions  for  the  cross-sectional 
area  computational  points.  The  boundary  Is  taken  to  be  at  the 
compressor  outlet  and  in  front  of  the  inlet  airfoil.  The  boundary 
conditions  are:  in  front  of  the  inlet  airfoil,  (dr/dz)a^onp.  stream  line 
=  0,  at  the  outlet  of  the  compressor,  (dr/dz)along  stream  nne«  0.  * 

3.  Computation  of  the  reference  attack  angle  is  based  on  the 
principle  in  Reference  [4],  but  the  original  airfoil  type  of  the 
sample  compressor  is  A40.  Due  to  the  lack  of  experimental  data, 
we  assume  temporarily  the  reference  attack  angle  to  be  2°  smaller 
than  that  of  the  NACA-65  series  of  airfoils. 


4.  The  angle  of  lag  the  reference  attack  angle  is  also 
computed  according  to  the  principle  in  Reference  [4],  When  the 
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attack  angle  deviates  from  the  reference  attack  angle,  the  change 
in  the  angle  of  lag  may  be  corrected  according  to  the  principles 
of  Reference  [5].  The  corrected  curve  is  as  shown  in  Figure  1. 

5.  Below  the  reference  attack  angle,  the  loss  calculation  is 
based  on  the  principles  of  Reference  [4]  —  namely,  calculating  the 
loss  with  the  pressure  expansion  factor  and  relative  blade  height 
as  parameters.  When  the  attack  angle  deviates  from  the  reference 
attack  angle,  the  loss  is  also  calculated  according  to  the  principles 
of  Reference  [5].  The  increase  in  entropy  may  be  calculated  according 
to  the  following  equation  from  the  total  pressure  loss  coefficient 

Ui  : 


Where  <*>  is  the  angular  velocity  of  the  rotor,  >•*  is 

the  relative  stagnant  temperature  of  the  preceding  cross-sectional 
area.  For  a  stationary  airfoil:  u  *  0,  we  get 


R 


ta{,  -  a  [l  -(l  +A— 


6.  The  effect  of  the  inner  and  outer  wall  boundary  layers  on 
the  flow  is  taken  into  account  by  coefficient  KB.  And  in  this 
calculation/  the  first  stage  KB  is  taken  to  be  0.99,  the  second 
stage  0.98,  for  the  3rd  stage  and  beyond  0.96,  and  kept  constant. 


7.  When  calculating  the  characteristics  properties  of  a  level, 
the  mean  value  of  the  level  is  calculated  according  to  the 
flow  rate  mean  value  for  parameters  at  various  streamlines. 


4.  COMPUTATIONAL  PROCESS 


For  our  sample  calculation,  since  (k  —  Data  9*  •  —  M‘m)  i 

far  less  than  1,  by  assuming  fit /dr  and  dt/  Qm  as  quantities  of  -.'r. 
same  order,  then  the  dt/dm  term  may  be  neglected.  After  Eq.  (6) 


40a 


is  dispersed  we  get : 


W\.i+X  -  WW  -  (2 gt/.li)T,Oi*t  -  *d 


where 


When  calculating  along  the  radial  direction,  we  first  assume 
the  axial  Velocity  for  the  root  cross-section  area,  and  gradually 
calculate  toward  the  head  section  to  find  Wz  and  p  of  each  point. 
When  the  whole  cross-section  is  calculated,  the  flow  rate  is 
obtained  by  integration.  If  the  flow  rate  does  not  satisfy  the 
requirement,  then  the  root  velocity  is  corrected.  The  correction 
formula  is 


wW  -  IT/AU  -  fcUGW  -  ct)/c,]> 

where  kQ  xs  the  iteration  factor  of  the  flow  rate  iteration.  When 
the  flow  rates  of  each  cross-sectional  area  are  all  satisfactory, 
then  stream  line  iteration  is  carried  out.  During  stream  line 
iteration,  the  correction  formula  for  the  stream  line  radial  position 
is 


rft°  -  -  'ft) 


Where  is  the  radial  position  of  the  stream  line  at 

the  t  Iteration,  W*:  the  current  value  of  the  stream  line  radial 

position,  ^7°:  the  radial  position  of  the  stream  line  at  the  t  +  1 

iteration,  and  k  is  the  stream  line  iteration  factor. 

*  r 


In  this  calculation,  we  take  the  relative  error  satisfied 
by  the  flow  rate  and  stream  line  to  be  10  . 

5.  DISCUSSIONS 

We  carried  out  computations  for  subsonic  compressor  stages  with 
the  computer  program.  Two  of  the  sample  calculations  have  the 
flow  through  part  shown  in  Figure  2  and  Figure  8.  The  low  speed 
experi  mental  results  for  the  stage  as  shown  in  Figure  2  are  given 
in  Reference  [6].  We  discuss  below  several  problems  based  on  the 
computed  results. 

1.  Loss  model:  To  obtain  results  in  conformity  with  practical 
results,  one  must  have  a  reasonable  loss  model.  Here  we  choose  the 
loss  model  given  in  Reference  [5].  The  loss  given  in  Reference  [7] 
and  [4]  is  larger  than  that  in  Reference  [5],  especially  in  the  head 
section.  For  single  stage  test  compressor  (channel  shown  in  Figure 
2),  the  stage  characteristics  calculated  according  to  the  loss  model 
in  Reference  [5]  and  [7]  are  shown  in  Figure  3.  Generally  speaking, 
under  an  ordinary  attack  angle,  although  the  effect  of  these  two 
loss  models  on  efficiency  and  pressure  ratio  is  not  large,  it  is 
still  apparent.  From  the  comparison  before,  we  know  that  the 
effect  of  the  loss  model  on  the  characteristics  is  weaker  than  that 
of  other  factors. 


Figure  1.  Calculation  of  the  angle  of  lag  under  varying 
working  conditions.  The  solid  line  is  from  Reference  [5]. 
The  dashed  line  is  the  correction  to  Reference  [5] 

m -0.5-u i  ,  the  subscript  r  is  the  reference  work 

condition  parameter. 
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Figure  2.  Diagram  of  single  stage 
test  compressor  flow  channel  and 
measured  cross-section. 

Key:  1 — measured  cross-section. 


2 .  The  reference  angle  of  attack  and  the  angle  of  lag  under 
the  reference  angle  of  attack:  The  effect  of  these  2  factors  on 
the  computed  results  of  the  stage  characteristics  is  shown  in  Figure 
3.  From  the  diagram,  it  can  be  seen  that  the  unsuitable  choice 

of  reference  angle  of  attack  has  a  greater  effect  on  efficiency, 
and  the  angle  of  lag  affects  both  pressure  ratio  and  efficiency, 
significantly,  and  its  effect  on  the  stage  characteristics  is  far 
greater  than  that  on  efficiency.  Hence,  only  by  determining  accurately 
the  reference  attack  angle  of  the  cascade  and  its  corresponding  angle 
of  lag  will  realistic  characteristics  be  found. 

3.  Convergence  of  stream  line  iteration:  For  different 
stream  line  iteration  factors,  the  number  of  iterations  satisfying 
the  same  accuracy  vary  widely.  For  simple  stage  test  compressor, 
when  \r  —  O.f  f  54  iterations  are  needed  to  satisfy  an  accuracy  of 
lO"24  while  for  kr  ■  0.4,  only  13  iterations  are  needed.  For  sub¬ 
sonic  stages  with  a  flow  channel  which  is  not  too  steep,  within 

a  wide  range,  all  iteration  factors  are  convergent,  but  when  the 


flow  channel  shape  varies  more  greatly,  the  Iteration  factor 
needs  to  be  chosen  small  to  guarantee  stability.  Thus,  selection 
of  the  stream  line  iteration  factor  should  be  determined  by  the  change 
in  stream  line  and  the  combination  of  various  parameters  (such  as 
slope,  curvature,  etc.).  Hence,  varying  the  iteration  factor  to 
affect  iteration  seems  to  be  more  effective. 
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Figure  3.  Effect  of  reference  attack  angle,  reference  lag 
angle,  and  selection  of  loss  model  on  stage  characteristics. 
Single  stage  test  compressor,  n=19000  revolution/min. 


4.  Effect  of  different  computational  methods  of  the  curvature 
and  slope  of  stream  lines  on  the  characteristics:  Besides  represent¬ 
ing  the  stream  line  with  a  cubic  representation,  it  can  also  be 
represented  as  a  multi-nomial.  If  we  still  use  the  2  cross-sectional 
areas  at  the  inlet  of  the  compressor  and  in  front  of  the  inlet  air¬ 
foil  as  the  reference  cross-sectional  areas,  then  in  the  case  of  our 
sample  calculation,  the  characteristics  obtained  by  calculating 
the  slope  and  curvature  of  the  stream  line  with  a  multi-nomial  yield 
the  same  result  as  with  a  cubic  representation.  The  difference  is 
negligible. 


5.  Effect  of  the  radial  distribution  of  curvature,  slope  and 
entropy  on  the  radial  equation  of  motion;  Figure  4  is  the  calculated 
results  of  the  radial  distribution  of  curvature,  slope,  and  entropy 
on  the  axial  velocity  distribution.  From  this  one  can  see  that 
generally  for  the  first  few  stages  of  the  subsonic  compressor, 
not  only  is  there  not  much  effect  on  the  axial  velocity  distribution 
due  to  the  non-uniformity  of  entropy,  but  for  ordinary  channels, 
there  is  not  much  effect  due  to  stream  line  curvature  or  slope 
either,  (generally  less  than  5%)- 
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Figure  4.  In  a  completely  radial  equation  of  motion,  the  effect 
of  the  variation  of  entropy,  curvature  and  slope  in  axial 
velocity  calculation,  (sample  computation  2) 

Key:  1 — radial  point;  2 — stationary  airfoil;  3 — moving  airfoil. 


In  addition,  if  only  the  axial  gap  dimensions  of  the  channel 
are  changed  and  not  the  radial  dimension  of  the  various  calculated 
cross-sectional  areas  or  the  geometric  parameters  of  the  airfoil, 
calculation  shows  that  there  is  very  little  difference  between  the 
characteristics  obtained  and  the  original  characteristics.  This 
also  shows  that  the  curvature  and  slope  of  the  stream  line  has  little 
effect  on  the  subsonic  stage.  Hence,  in  simulation  testing,  to 
make  the  measurement  more  convenient,  increasing  the  axial  dimensions 
will  still  give  usable  experimental  results  for  the  subsonic  stage. 
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Figure  5.  Comparison  of  perfect  radial  equation  of  motion 
and  simple  radial  equation  of  motion  —  with  the  outlet 
entropy,  axial  velocity  and  change  in  attack  angle  of  sample 
calculation  2. 

Key:  1 —  — o— o—  simple  radial  equation  of  motion 

—A— A—  perfect  radial  equation  of  motion; 

2 — stationary  airfoil;  3 — radial  point. 


6.  Comparison  of  simple  radial  equation  of  motion  and  perfect 
radial  equation  of  motion:  For  sample  calculation  2,  Figure  5  shows 
the  results  calculated  for  the  simple  radial  equation  of  motion  and 
perfect  radial  equation  of  motion,  including  the  entropy  behind 
a  moving  airfoil  and  behind  a  stationary  airfoil,  axial  velocity 
and  the  attack  angle  of  a  moving  airfoil  and  of  a  stationary  airfoil, 
etc.  From  the  diagram  ,  it  can  be  seen  that  there  is  definitely 
a  difference  in  the  2  attack  angles  calculated.  This  is  also  reflected 
in  the  radial  distribution  of  the  axial  velocity  and  the  stage 
characteristics.  But,  if  the  slope  of  the  passage  way  is  not  large, 
this  difference  is  small.  But,  if  the  slope  of  the  passage  way 
is  large,  then  it  is  better  to  calculate  according  to  the  perfect 
radial  equation  of  motion.  Figure  6  gives  the  difference  between 
the  characteristics  for  these  2  methods  of  calculations. 


1.26 

1.24 


Figure  6.  Comparison  of 
characteristics  for  the 
perfect  radial  equation  of 
motion  and  the  simple  radial 
equation  of  motion. 

(sample  calculation  2  — 
n  =  8600  rev/min) 

Key:l  -a-a-  perfect  radial 
equation  of  motion, 

2  -o-o-  simple  radial 
equation  of  motion. 


6.  COMPARISON  WITH  TEST  RESULTS 
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Figures  6-8  gives  the  diagrams  of  the  stage  characteristics 
and  stream  line  in  the  2  sample  calculations.  Figure  9  gives  the 
air  stream  exit  angle  inlet  leading  airfoil,  moving  airfoil  and 
stationary  airfoil.  Figure  7  and  Figure  9  also  show  the  experimental 
results.  From  the  diagrams,  it  can  be  seen  that:  in  the  vicinity 
of  the  reference  attack  angle,  the  calculation  agrees  with  the 
experiment,  but  when  we  deviate  from  the  reference  attack  angle 
region,  there  is  a  larger  difference  between  calculation  and 
experiment  in  the  characteristic,  especially  in  the  pressure  ratio. 
For  positive  attack  angle  ,  the  calculated  pressure  ratio  is  higher 
than  the  experimental  value.  One  may  regard  this  as  due  to  the 
angle  of  lag  and  loss  under  changed  working  conditions  in  Reference 
[5],  especially  due  to  the  unreasonable  calculation  of  the  angle 
of  lag.  We  made  corrections  to  the  curve  of  Reference  [5]  based 
on  this  (see  Figure  1)  and  increased  the  change  in  the  angle  of 
lag  below  the  positive  attack  angle.  The  characteristics  calculated 
from  the  curve  of  corrected  lag  angle  changes  are  also  plotted  in 
Figure  7  which  shows  good  agreement  with  experiment.  From  Figure  9, 
we  see  that  in  the  vicinity  of  the  design  point,  the  calculated 
lag  angle  value  also  agrees  with  the  experimental  value. 


c\T\  p* .  io’-  a\T\  v* .  10- 

n=100o0|f/^j*  n=ll000H.  _/)■ 

Figure  7.  Characteristics  of  single  stage  test  compressor 
(n =10000  and  n=l4000  rev/min. ) 

-  perfect  radial  equation  of  motion. 

-  simple  radial  equation  of  motion. 

Key:  1 — •  calculated  point ;  2 — x  calculated  point,  from 
corrected  curved  in  Figure  1;  3 — o  experiment  point; 

4 — •  calculated  point;  5 — xcalculated  point,  based  on 
corrected  curve  in  Figure  1;  6 — o  experiment  point. 


0  0.1  0.2 
i(m) 


Figure  8.  Channel  shape  and  stream  line  distribution  in 
sample  calculation  2.  (n=8600,  G=48.5  kg/sec) 

- o -  experimental  *>  -  0.781 

- • -  calculated  G  =  5.0  kg/sec. 
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ir  !U‘  2lT  10*  20*  30"  40*  50‘ W  10'  20‘ 


Figure  9.  Single  stage  test  compressor  air  stream  exit  angle  — 
comparison  between  calculated  and  experimental  results. 

>'r.*10000  rev /min.  ) 

" .  CONCLUSIONS 

1.  For  subsonic  stage  with  small  channel  slope,  there  is  not 
much  difference  in  the  calculated  results  between  a  perfect  radial 
equation  of  motion  and  simple  radial  equation  of  motion.  But  when 
the  slope  is  larger,  one  should  calculate  with  the  perfect  radial 
equation  of  motion. 

2.  For  the  calculation  of  compressor  characteristics, 
the  degree  of  accuracy  of  the  reference  angle  of  attack  and  its 
corresponding  angle  of  lag  is  more  important  than  the  choice  of 
loss  model. 

3.  A  problem  requiring  immediate  solution  in  characteristic 
calculations  is  the  method  of  calculating  the  angle  of  lag  and 
loss  under  changing  working  conditions. 

4.  In  order  to  apply  3  dimensional  flow  theory  effectively 
to  engineering,  the  most  pressing  task  currently  is  to  strenthen 
experimental  research.  It  will  not  only  provide  experimental 
data  for  the  proposed  physical  models,  but  also  will  give  the 
criteria  for  judging  whether  a  physical  model  is  good  or  tad. 
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APPLICATION  OF  THE  FINITE  ELEMENT  METHOD  TO  THE 
SOLUTION  OF  TRANSIENT  TWO-DIMENSIONAL 
TEMPERATURE  FIELD  FOR  AIR-COOLED 
TURBINE  BLADE 

Guo  Kuan-liang  Ge  Xin-shi  Sun  Ziao-lan 
(The  Chinese  University  of  Science  of  Technology ) 


Abstract 

To  assure  air-cooled  turbine  blade  of  being  able  to  operate  safely  and  reliablely 
during  continuous  variation  of  mechanical,  particularly  thermal  loads,  it  is  necessary  to 
calculate  the  transient  temperature  field  of  blades.  In  this  paper  the  finite  element 
method  for  computation  of  transient  two-dimensional  temperature  distribution  is  dis¬ 
cussed.  the  calculation  of  transient  temperature  field  for  a  middle  cross  section  of  ;iir- 
cooled  blade  is  made  and  discussed. 


1.  INTRODUCTION 

The  adoption  of  air-cooled  blades  has  Increased  the  inlet 
temperature  of  turbines  and  provided  a  powerful  measure  in  improving 
the  performance  of  jet  engines.  To  make  sure  that  the  air-cooled 
blades  work  safely  and  reliably  at  high  temperatures,  high  pressure, 
burning  gas,  the  calculation  of  the  temperature  field  is  a  vital 
link  in  the  research  and  development  of  air-cooled  blades.  Especially 
when  the  engine  is  working  under  acceleration  or  deceleration,  the 
mechanical  and  thermal  loads  of  the  blades  change  drastically  and 
the  blades  are  working  under  extremely  severe  conditions.  Hence, 
the  calculation  of  a  transient  temperature  field  is  even  more 
important . 

* 

This  paper  has  been  presented  at  the  Second  National  Engineering 
Thermal  Physics  Conference  of  Hangzhow  in  November,  1978. 
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As  far  as  we  know,  there  is  a  large  quantity  of  work  on  the 
finite  element  method  for  solving  problems  involving  thermal  conduc¬ 
tion  of  solids  in  foreign  countries,  but  very" little  has  been 
reported  on  using  the  finite  element  method  to  calculate  the  trans¬ 
ient  temperature  field  of  air-cooled  blades  and  the  treatment  of 
boundary  parameters. 

In  recent  years,  works  have  been  started  domestically  in  apply¬ 
ing  the  finite  element  method  to  calculate  the  steady  state  tempera¬ 
ture  field  of  air-cooled  blades.  In  Reference  [1],  a  method  of 
air-cooled  blade  boundary  parameter  treatment  is  presented,  and  a 
special  computer  program  in  ALGOL-60  language  on  the  DJS-21  computer 
is  given.  Reference  [2]  summarized  how  to  establish  the  variational 
principles  thorugh  analyzing  the  physical  phenomena  from  which  the 
equations  for  the  finite  element  method  of  solving  the  steady  state 
and  transient  thermal  conductivity  problems  are  derived.  Reference 
[3]  proposed  the  computational  method  for  the  transient  temperature 
field  of  the  air  stream  in  the  passageway  of  air-cooled  blades. 

In  this  paper,  we  have  compared  the  effects  of  various  time 
step  sizes  on  computational  results  and  presented  a  FORTRAN  program 
on  the  DJS-8  computer  for  calculating  the  2-dimensional  transient 
temperature  field  of  air-cooled  blades. 


2.  FINITE  ELEMENT  METHOD  FOR  2-DIMENSIONAL  TRANSIENT  THERMAL 
CONDUCTIVITY  PROBLEMS 


For  2-dimensional  transient  thermal  conductivity  problems  with 
boundary  conditions  of  the  3rd  kind,  the  differential  equation  is 


0 


on  region  R 


(1) 


The  boundary  condition  is 


/l4l 


dT/dn  -  (a,/l)(T,  -  T)  on  boundary  r  (2 ) 

For  air-cooled  blades,  both  the  inner  and  outer  walls  will 
exchange  heat  convect ionally  with  the  fluid.  The  corresponding 
boundary  condition  may  be  written  as 

dT  fO,/A)(T.,  -  T)  on  the  burning  gas  side 
d«"  ™  Ko^/iXT,  —  T)  on  the  cool  gas  side 


The  initial  condition  is 


where:  p*e>l  are  respectively  the  density,  specific  heat  and 
average  thermal  conductivity  of  the  blade  material.  <*»»<*. 

are,  respectively,  the  convectional  heat  exchange  coefficients 
of  the  burning  gas  and  the  cool  gas  on  the  blade  wall.  T.ttTm 
are,  respectively,  the  absolute  thermal  recovery  temperatures  of 
the  burning  and  cooled  gas.  Since  the  stream  velocity  of  the  cool 
gas  is  rather  low,  therefore  where  the  subscript  s  denotes 

air  streams. 

To  use  the  finite  element  method  for  solving  the  transient 
temperature  field,  we  must  first  obtain  the  equivalent  starting 
equations  of  the  finite  element  method  for  the  above  problem.  This 
is  usually  done  through  the  mathematical  treatment  of  the  relatively 
simple  method  of  the  weighted  residual  (e.g.,  Galerkin  method) 
on  from  the  variational  principle  which  is  more  intuitive  in 
its  physical  meaning  [2], 


Here,  we  present  the  result  derived  with  the  Galerkin  method 


[4,5] 


dx  9x 

+  (  8T*-(.T- 
Jr  3L 


(4) 
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1.  Spatial  discretization: 

When  we  discretize  equation  (4),  we  assume  that  the  thermal 
expansion  of  the  solid  may  be  neglected. 


We  divide  the  solution  region  0  into  NE  triangular  surface 
elements,  and  represent  respectively,  with  NP  and  NL,  the  number 
of  nodal  points  (i.e.,  the  total  number  of  vertices  of  the  triangles) 
and  the  line  element  on  the  boundary  r  .  Consider  the  linear 
insertion  function  of  an  arbitrary  triangular  surface  element 
(with  vertices  i,j,m) 

T  —  (5) 


where  Nt,  and  Nm  are  the  shape  function  and 

2A  {  A v»)/2  (6) 

«|  —  *ijm  —  b,  —  y,  —  ym;  <■(  —  *.  —  x, 

•i  -  *.y,  -  b,-ym-y,i  e,  — x,  -  x„  (7) 

a»  ■“  *i7i  —  *iyii  bm  —  JP|  —  yr,  fm  —  X,  —  X, 


Take  *'T  with  the  same  functional  form  as  T,  i.e. 

V/ 


(8) 


where  are  respectively,  the  value  taken  by  the  arbitrary 

function  8T  at  nodal  points  i,j,  and  m.  Similarly  for  a  line 
element  on  the  boundary  r  (let  its  end  points  be  i  and  j)  we  may 
take 


(9) 


/142 


By  using  eq.  (5)-(9),  we  may  discretize  the  terms  in  eq.  (4)  as 


IJjrf  r'"'-TCI' 

j r8T  ?  (T  ~  T^dl  “  2  *  j(T  -  7,)J/  -  f(//r-  *) 


and  also  get 


Ct  +  KT-K-o  (10) 

where  #C  —  1C*  +  //. 

After  eq.  (4)  Is  spatially  discretized,  we  obtain  a  set  of 
equations  including  NP  first  order  ordinary  differential  equations 
as  mentioned  above.  The  appearance  of  the  matrix  C  is  related  to 
the  thermal  capacity  term;  the  in  the  K  matrix  has  nothing  to 
do  with  time  while  H  is  related  both  to  space  and  to  time  (it 
contains  the  local  thermal  radiation  coefficient  and  thermal 
conductivity);  R  is  a  column  vector,  and  is  homogeneous;  the 
column  vector  T  is  the  temperature  function  yet  to  be  found. 

2.  Time  discretization 

To  convert  (10)  into  linear  algebraic  equations,  we  may  use 
different  methods  of  time  discretization.  Three  discretization 
schemes  are  presented  here: 

(1)  By  using  the  method  of  a  weighted  residual  with  linear 
approximation  for  the  temperature  function,  etc.,  it  is  easy  to 
show  that 


(O/Aixc;  2Ci)  +  x,  +  jir,ir, 

-  [(2/aocc;  +  2Ci)  -  k,  -  jr,jr,  +  o*  +  4*) 


where  — *»*-*•  is  the  time  step. 


(2)  By  using  the  method  of  weighted  residual  with  quadratic 
approximation  for  the  temperature  function,  etc.  [6,7],  we  can 
derive 


MuljT.l  r,l 

~Mu  A/JlrJ  lC|  —  Ml,  7*»J 


(12) 


where 


M  it  “  KLiiKl  4*  K/40hH  +  K  A I  t!H%  +  KA2hH\  +  C EaC/At/2 
C,  -  -(££.,*,  +  RE, +  REM 


KL  - 

'224  28  281 

■  28  56  -14-1 

KAO  — 

r  .« 

-8 

1 

L— 8 

-3 

KA\  - 

192  16  161 

■  16  20  — 8J 

KA2  - 

f“ 

20 

(.20 

39 

C£-  I 

0  140  —1401 
-140  105  35J 

RE  - 

f  " 

224 

1 

L— 14 

28 

(3)  C-N  scheme  [8,9,10] 

Let  us  assume  that  within  the  time  step  &fc ,  the  derivative 
T  of  temperature  wrt  time  varies  linearly  with  time,  then 

(7*.  -  T,)/£u  -  (f.  +  tx)/l 

Assume  also  that  (pf/AV(<w/i>,  —  i,  i.e.,  —  I  (one  may 

not  make  this  assumption)  the  subscript  1  means  the  inverse  matrix 
of  the  labeled  matrix.  Then  we  may  transform  eq.  (10)  into 


1(2 C,/A<)  +  JC.ir,  -  [(2C1/A1)  -  K,]T%  +  R.  +  *, 


(13) 


It  is  worth  while  to  point  out  that  since  in  eq.  (10),  C,K,  and 
R  are  all  functions  of  temperature,  it  is  a  set  of  non-linear 
ordinary  differential  equations.  When  making  use  of  the  C-N 
scheme,  we  have  considered  the  non-linearity  of  the  equations. 
The  result  obtained  (13)  is  an  improvement  over  the  usual  linear 
approximations . 
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(4)  Comparison  of  the  3  schemes  above 

The  results  of  3  discretization  processes,  i.e.,  eq.  (11), 

(12),  and  (13)  are  all  implicit  and  are  unconditionally  stable. 

When  the  3  schemes  are  composed,  we  may  for  convenience 'assume 

eq.  (10)  is  linear  and  homogeneous.  Then  it  is  not  difficult 

to  prove  that:  Eq.  (11)  is  accurate  to  order  0(h)  (h  being  the  time 

step  length)  and  eq.  (12)  and  (13)  are  both  accurate  to  order 
2 

0(h  )  (better  than  eq.  (11),  but  since  the  latter  is  much  simpler 
to  calculate  than  the  former,  we  shall  use  eq.  (13)  in  our  calcula¬ 
tion. 


3.  COMPUTATION  OF  BOUNDARY  PARAMETERS 

In  order  to  calculate  the  various  instantaneous  C,  K,  and  R 
in  eq.  (13),  and  then  solve  for  the  transient  temperature  field  of 
the  air-cooled  blade,  we  must  first  know  the  instantaneous  values 
of  the  boundary  parameters  a«,  <*.»  T4t  and  and  the  average 

thermal  conductivity  X  and  specific  heat  C  of  the  blade  material. 

In  this  work,  we  calculated  a  transient  process,  starting  with  the 
initial  state  of  the  air-cooled  blade  (assumed  to  be  a  steady  state) 
and  gradually  accelerating  for  a  given  period  time  until  the 
parameters  of  the  burning  gas  and  cool  gas  stream  no  longer  change 
and  finally  after  a  period  of  "stagnation"  reaching  the  steady  state 
at  the  end.  Of  course,  the  finite  element  method  and  the  method  of 
treating  boundary  conditions  as  mentioned  above  are  also  applicable 
to  other  state-changing  situations.  In  the  following,  we  shall 
briefly  introduce  the  method  of  treating  the  parameters.  For  a 
detailed  calculation,  see  Reference  [3]. 

1.  If  we  denote  by  ■  t,  (the  subscript  i  indicates  the  parameter 
under  discussion)  the  ratio  of  the  various  instantaneous  parameters 
of  the  burning  gas  and  cool  gas  stream  to  the  corresponding  parameter 
at  the  final  steady  state,  then  and  can  be 

conveniently  obtained  from  the  known  value  of  5-  and  the  steady 
state  values  of  and  r»»«  already  obtained  in  Reference  [1], 
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and  there  is  no  need  to  go  through  the  complicated  computation 
according  to  the  original  standard  formulas  one  by  one. 

For  example: 

«.  -  •  SimCSwt/SH)w  " 

<*, — «« • 

2.  The  thermal  capacity  of  the  blade  must  be  taken  into  account 
in  calculating  the  transient  temperature  cf  the  cool  gas  stream  along 
the  blade  wedge  direction.  Here  the  basic  equations  treating  by 
sections  according  to  one  dimensional  problem  should  be 

4 

—  ?i)  “  T*)  4*  7200rfK(T|  —  i*»)/ A* 

*«Ptk(T.  -  ?„)  -  G'tC^TU  -rj  *  -  1  ,2,  3,  4 
f.k  -  CTU  +  T%)/2  -t- 1,2,  3,4 

*-4,3,2 

where:  Tck  is  the  cool  gas  stream  temperature  in  the  kth  channel; 
y  is  the  specific  gravity  of  the  blade  material;  V  is  the  volume 
of  each  blade  section;  P„  f,k  are  respectively,  the  outer  area 
of  the  blade  and  the  inner  area  of  the  kth  channel;  f„  is 
the  average  temperature  of  each  blade  section  at  the  beginning 
instant  of  the  time  step  At,  its  initial  value  being  calculated 
from  the  initial  steady  state;  the  subscripts  i  and  o  denote, 
respectively,  the  parameters  at  the  outlet  and  inlet  of  each 
section. 

Since  fm  and  T«  of  the  above  equations  are  already  known, 
by  solving  the  above  equations,  we  can  then  get  the  temperature 
variation  of  the  cooling  gas  stream  in  the  passages  of  the  blade 
at  each  instant,  and  also  the  average  temperature  of  the  /iM 

blade  cross-section  from  which  the  corresponding  -p%  *  and  i, 
may  be  found  and  and  computed. 


4.  THE  COMPUTER  PROGRAM  AND  ITS  EXPLANATION 


From  eq.  (13'  it  is  not  difficult  to  see  that  although  TQ.is 
given  by  the  initial  conditions  (3)  and  the  initial  state  Rq  and 
Kq  can  be  calculated  from  known  parameters,  yet  R^ ,  and 
are  both  functions  of  T^.  Hence,  it  is  necessary  to  use  the 
iteration  method  to  solve  for  T^.  Since  we  have  already  obtained 
the  cross-sectional  average  temperature  of  the  blade  at  each  instant 
when  we  calculate  Tck,  it  may  be  used  as  a  very  good  approximation 
of  the  cross-section  average  of  T^.  With  it  we  can  calculate  a 
prior  R^,  and  and  the  R^,  and  at  each  instant.  Eq.  (13) 
may  then  be  transformed  to  the  following  form 


KT-R 

and  may  be  solved  by  using  the  direct  matrix  analysis  method.  This 
will  not  only  satisfy  the  demand  of  engineering  computations,  but 
also  reduce  greatly  the  amount  of  computation. 

Because  of  the  limitation  of  space,  we  shall  omit  the  computational 
flowchart,  the  FORTRAN  program  and  explanations  in  this  paper. 


5.  COMPUTATIONAL  RESULT  AND  ANALYSIS 

In  this  work,  we  have  calculated  the  temperature  distribution 
of  the  blade  at  each  instant  of  an  accelerated  working  state.  For 
an  intuitive  understanding.  Figure  1  and  2  show  the  instantaneous 
temperature  distribution  at  5  different  instants  at  the  periphery 
on  the  burning  gas  side  of  the  blade  middle  cross-section  (i.e.  , 
blade  back  and  blade  bowl).  In  Figure  3,  the  time  variation  of 
the  average  temperature  of  the  middle  cross-section  calculated 
by  the  finite  element  method  is  plotted.  Also  for  comparison, 
the  cross-sectional  average  temperature  curve  obtained  from  one 
dimensional  sectional  treatment  is  also  plotted. 


1.  The  point  of  maximum  temperature  on  the  air-cooled  blade 
is  located  at  the  boundary  of  the  tail  section  on  the  burning  gas 
side,  and  the  point  of  lowest  temperature  is  located  at  the 
boundary  of  the  internal  cooling  orifice.  This  is  quite  reasonable. 
The  isothermal  distribution  curve  at  each  instant  also  agrees  with 
physical  law. 

It  can  be  seen  from  Figure  1  and  2  that  the  temperature  at  the 
front  edge  and  tail  edge  is  higher  and  that  at  the  middle  of  the 
blade  back  and  bowl  is  relatively  lower.  When  the  blades  are 
continuously  accelerated,  their  temperatures  continue  to  rise.  The 
tendency  of  the  temperature  distribution  at  different  instants  is 
basically  the  same. 

2.  From  Figure  3,  comparing  the  cross-sectional  average 
temperature  calculated  from  the  finite  element  method  and  that  from 
one  dimensional  treatment  by  sections,  the  maximum  difference  is 
below  2%.  This  shows  that  when  solving  eq.  (13),  by  calculating 

a  priori  C1,  K1  and  with  the  various  instantaneous  cross- 
sectional  average  temperatures  obtained  by  sectional  treatment  and 
reducing  the  non-linear  problem  to  linear  equations,  it  is  definitely 
possiole  to  provide  the  accuracy  required. 

3.  To  compare  the  effect  of  various  time  step  length  At,  in 

this  calculation,  we  use  respectively,  0.02,  0.5,  1  and  2  seconds  as 
the  time  step  length.  The  results  indicate  that  the  time  step  lengths 
used  have  little  effect  on  the  cross-sectional  average  temperature.  For 
example,  when  A/  — 0.02  sec,  1008.2  £  at  16  seconds  which  for 

A<  —  2  second,  T»—  1011.4  K  at  16  seconds.  Thus  when  the 
time  step  length  is  increased  100  times,  the  blade  cross-sectional 
average  temperature  only  varies  by  3.2  degrees.  The  relative  differ¬ 
ence  is  only  0.3%  (the  other  cases  are  even  better  than  this). 

4.  Similar  to  the  situation  of  finding  the  2  dimensional 
steady  state  temperature  field  of  air-cooled  blades  with  finite 
element  method  [1],  the  accuracy  of  the  transient  temperature  field 


calculation  is  also  determined  basically  by  the  calculation  of  the 
boundary  parameters.  Hence,  accurate  determination  of  the  local 
convective  heat  exchange  coefficient  «f  on  the  burning  gas  side, 
the  convective  heat  exchange  coefficient  a,  of  each  internal 
cooling  passages  and  the  cool  air  flow  rate  c,  of  the  air-cooled 
blade  requires  further  investigation. 

In  completing  this  work,  we  have  received  the  earnest  support 
of  the  Research  Institute  of  Shenyang  Aeronautic  Engine  Company. 
Comrades  Li  Yishen,  and  Huang  Kecheng  of  our  university  also  gave 
us  helpful  support.  We  would  like  to  express  our  gratitude  to 
all  of  them. 
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Abstract 

This  article  discussed  t li«*  physical  imulel  for  thermal  diffusivity  measurement  by 
pulsed  laser  technique,  and  described  the  fundamentals  of  data  processing  and  contr.d 
by  a  computer.  Several  boundary  conditions  must  be  satisfied  in  the  measurement. 

The  thermal  diffusivities  of  standard  specimens - Armco  iron  and  a-Al.O,,  as  well 

as  other  materials  were  determined  by  the  apparatus.  Data  obtained  from  the  stan¬ 
dard  specimens  were  compared  with  recommended  values  in  the  literature  and  found 
to  be  in  good  agreement.  The  r.  m.  s.  errors  of  Armco  iron  and  a-Al,0,  specimen 
measurements  amount  to  -3.2—  +6.3%  and  -2.4  — +5.6%  respectively.  This  ap¬ 
paratus  features  automatic  and  lngh  speed,  measurement,  and  is  ricommeuded  to  be 
operated  in  the  temperature  range  of  30<>  to  1800°<\  The  interval  between  laser  erais- 
>i<  i  and  display  of  result  is  approximate! .  5  seconds. 


1.  INTRODUCTION 

The  pulsed  laser  technique  to  measure  thermal  diffusivity 
A  as  developed  in  the  sixties  has  overcome  the  limitations  of  other 
methods,  and  provides  the  advantages  of  small  sample  size,  high 
speed,  high  precision  and  is  applicable  to  the  measurement  of 
a  wide  variety  of  materials.  Hence,  it  attracts  cl.ose  attention 
and  develops  very  rapidly.  According  to  reports  [1,2],  75%  and 
more  of  the  thermal  diffusivity  data  have  been  measured  with 
this  method  after  1973.  In  the  seventies,  after  the  incorporation 
of  computer  technique,  the  precision  and  speed  of  laser  thermal 
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diffusivity  measurement  apparatus  have  pronouncedly  Improved  [2]. 

Our  Institute,  together  with  the  Shanghai  Laser  Station, 

Shen  Yang  Metal  Institute,  Academia  Sinica  and  Shanghai  Steel 
Institute,  jointly  researched  and  developed  the  laser  thermal 
diffusivity  measurement  apparatus  in  1973.  In  April  of  1978,  our 
institute  also  successfully  developed  a  computerized  laser  thermal 
diffusivity  measurement  apparatus  and  improved  its  various 
characteristics . 


2.  METAL  PHYSICS  MODEL  AND  PRINCIPLE  OF  MEASUREMENT 

The  physical  foundation  of  the  steady  state  method  and  transient 
method  for  measuring  the  thermal  diffusivity  X  is  the  differential 
equation  of  thermal  diffusion.  This  equation  may  be  solved  by 
using  various  specific  boundary  conditions.  The  solution  obtained 
is  generally  related  to  the  thermal  properties  of  matter  [4], 

The  physical  model  of  the  pulsed  laser  measurement  technique 
for  X  is  to  radiate  normally  on  the  front  surface  of  a  thin  circular 
disc  (the  sample)  insulated  at  the  periphery,  with  a  uniform  laser 
pulse  beam  and  then  by  measuring  the  temperature  rise  curve  on  the 
back  side  of  the  sample  under  one  dimensional  heat  flow  conditions, 
to  find  the  temperature  diffusivity  a  from  whicb  the  x  value  may 
be  deduced. 


Let  the  initial  temperature  distribution  at  any  point  X  on 
a  circular  disc  with  sample  thickness  L  be  T(x,0),  then  the  temperature 
distribution  T(x,t)  [51  at  any  time  t  Is 


«».«>  -  \  j><„  o)  *  S>. «) 


nwx  j 
co* -  dm  . 
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(1) 


At  the  instant  /*(<•-*  0),  when  the  sample  absorbs  the  laser  /14-8 

the  temperature  distribution  at  X  within  a  very  small  distance  E 
from  the  front  surface  (*  —  o)  of  the  sample  is 
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(2) 


—  Q/D  •  e  •  g,  13  0  <  x  <  f ; 

T’C**  <•.)  ~  0,  g  <  x  <  L. 

(3) 

where  D  and  c  are  respectively  the  density  and  specific  heat  of 
the  sample,  Q  the  radiation  intensity  of  the  laser  pulse  absorbed 
by  the  sample  and  T(x,to)  is  the  temperature  rise  of  the  sample 
relative  to  the  experimental  store  temperature. 


Substituting  eq.  (2),  (3)  into  (1),  because  r» °»  T(*>  *•) 33  °) 


then 


n“ 0  “  ztott + i  2»"( 
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—  nnx 
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—  sVW 
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-  -  *  — sin—* 

L  D  '  c  •  g  nn  L 
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(4) 


To  simplify  the  temperature  measurement,  we  choose  for  the  back 
side  of  the  sample  (x=L)  the  following  special  case: 


X  —  L, 


n  xx 


“*  CO*»«r  «-  (  “  1 )’; 


(5) 


X  — o,  da  !a^!3 
4.  *  L  L  * 


(6) 


Substituting  (5),  (6)  into  (4),  we  get 


+  2  §<->•  •  -(=*4^)1  • 


(7) 


To  simplify  matters,  define  2  dimensionaless  parameters  to  and 
V(L,t)  and  let 


<« 

V(4.,i) 


L>  * 

nL,i) 

r«  * 


(8) 

(9) 


where  T^  is  the  maximum  temperature  rise  after  the  sample  is 
radiated  by  the  laser,  and 


—  Q/L  '  D  '  c 


(10) 


Substituting  (8)— (10)  into  (7),  we  get 

V(L,  0  —  l  +  2  2  (*" l)*exp(— *2u).  (11) 


with  V(L,t)  and  u  as  coordinate,  then  the  curve  in  Figure  1 
represents  eq.  (11).  From  eq.  (11)  and  Figure  1,  we  can  fix  the 
value  of  a.  When  ,  w—1.38  ,  then  Reference 

[8]  becomes 


_  0.138  •  L 1 

a  h  — ^  — 


(12) 


where  t?  is  the  time  required  for  the  back  side  temperature  of  the  sample 
to  reach  \  its  maximum  value,  and  may  be  found  from  the  temperature 
increase  curve  of  the  back  side  of  the  sample  from  which  thermal 
conductivity  a  may  be  found. 


The  boundary  conditions  required  by  eq.  (12)  must  be  satisfied 
during  the  measurement:  we  must  establish  a  one  dimensional  heat 
flow  along  the  sample  axis;  the  laser  pulse  time  must  be  far  smaller 
than  the  time  for  the  sample  to  reach  thermal  equilibrium;  the  laser 
pulse  should  be  uniformly  distributed;  the  heat  loss  of  the  same 
should  be  minimized.  Otherwise  corrections  must  be  made.  Under 
the  condition  of  second  apparatus  permit,  the  temperature  rise  on 
the  backside  of  the  sample  should  be  as  small  as  possible  so  that 
a  may  be  treated  as  constant. 

After  a  is  measured,  we  may  find  X  from  the  following  equation.  / 

i-a-fD  (13) 

The  error  of  C  value  as  measured  by  the  ice  calorimeter  or 
copper  calorimeter  is  generally  less  than  * — •  -5 %U*T’  ;  hence, 

according  to  the  usual  way,  the  value  of  C  may  be  found  usually 
in  handbooks. 
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Figure  1.  Curve  of  eq.  (11). 


3.  STRUCTURE  OF  THE  APPARATUS  AND  THE  COMPUTER  CONTROLLED 
MEASUREMENT  SYSTEM 

The  apparatus  basically  consists  of  4  parts  (see  Figure  2 
and  Figure  3)- 

1.  Laser  source:  a  Nd  glass  laser  with  xenon  lamp  excitation 
is  used.  The  wavelength  is  1.06  u,  Each  pulse  is  about  0.6  ms  with 
an  energy  of  about  20  Joules. 

2.  Vacuum  carbon  tube  oven:  the  oven  is  insulated  by  carbon 
fibre  and  heated  with  a  10  KVA  high  current  heater.  The  oven 
cavity  contains  graphite  sample  stand.  A  three-support  point  Zr02 
ring  is  used  in  the  oven  support  to  fix  the  sample  and  to  reduce 
heat  loss.  The  vacuum  chamber  is  water  cooled  and  is  capable 

of  automatic  elevation. 

3.  Opticoelectric  measuring  system.  The  sample  back  side 
temperature  rise  is  measured  with  a  light  sensitive  PbS  resistor. 

It  is  installed  in  a  blackened  tube  with  adjustable  grating. 
Infrared  lens  with  adjustable  focus,  and  filter.  The  optically 
sensitive  resistance  is  one  arm  of  a  bridge.  When  the  laser 
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Figure  2.  Picture  of  the  apparatus. 


Figure  3.  Block  diagram  of  apparatus. 

Key:  1 — parallel  lamp;  2 — laser;  3 — sample;  4 — infrared 
detector;  5 — electric  bridge;  6 — laser  power  supply;  7 — oven 
tube;  8 — vacuum  system;  9 — amplifier;  10 — oscilloscope; 

11 — CPU;  12 — ALU;  13 — A  to  D  converter;  14 — memory;  15 — 
quartz  clock;  16 — 38-10A  computer;  17 — output  device; 

18 — input  device;  19 — result  of  computation. 


radiates  on  the  front  side  of  the  sample,  the  temperature  rise 
on  the  back  side  increases  the  light  sensitive  resistance.  The 
bridge  becomes  unbalanced  and  outputs  a  small  voltage  on  the  load 
which  after  direct  current  differential  amplification  is  sent  to 
the  oscilloscope  with  the  time  marker.  The  temperature  increase 
curve  is  photographed  on  ultra-violet  film  or  the  signals  can  be 
sent  to  the  computer.  In  both  cases,  o  value  may  be  calculated. 


/150 
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4.  Computer  control  system:  we  used  a  JS-10A  computer  for 
the  control.  It  has  a  computational  speed  of  30000  computations/ 
second  and  a  memory  of  4k.  To  satisfy  the  special  requirements 
of  thermal  diffusivity  apparatus  control,  we  made  some  modifications 
to  the  computer  system  as  well  as  the  electronic  circuitry.  The  main 
control  program  designed  by  us  consists  of  5  sub  programs: 

1.  Recharging  and  starting  the  pulsed  laser  power  supply; 

2.  Collecting  and  storing  sample  back  side  temperature  rise 
curve  through  A  to  D  converter; 

3.  Processing  and  analyzing  the  temperature  rise  curve; 

4.  Calculating  the  values  of  t\  and  a; 

5.  Automatic  printout  of  experimental  results  and  all  the 
original  data. 

Before  testing,  the  above  programs  are  input  into  the  computer 
from  paper  tape;  then  the  data,  and  sample  thickness  are  entered 
into  the  computer  display.  The  main  switch  of  the  computer  is  then 
pressed  and  in  5  seconds  the  laser  is  charged  up  and  activated. 

In  another  5  seconds  or  so,  the  experimental  results  will  be 
displayed  on  the  printer. 


4.  THEORETICAL  ANALYSIS  OF  THE  VARIOUS  BOUNDARY  CONDITIONS  IN 
THE  TEST 

1.  Criteria  on  the  degree  of  satisfaction  of  boundary 
conditions. 

(1)  Analyzing  the  characteristics  of  the  temperature  rise  curve 
on  the  back  side  of  the  sample:  Figure  1  is  the  sample  temperature 
rise  curve  under  the  circumstance  that  the  experiment  satisfies 
all  the  boundary  conditions.  Compared  to  the  actually  measured 
temperature  rise  curve,  if  the  latter  reaches  the  final  equilibrium 
temperature  very  slowly  after  a  rapid  temperature  rise,  it  may  be 
due  to  the  non-uniformity  in  the  radiation  heating  of  the  sample  front 
side.  The  result  will  lead  to  a  higher  value  for  a.  If  there  is  a 
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serious  radial  heat  flow  in  the  sample,  then  the  temperature 
rise  curve  will  be  twisted.  If  the  temperature  rise  curve  decreased 
with  a  large  slope  immediately  after  reaching  the  maximum,  then  it 
indicates  a  large  radiation  heat  loss,  causing  a  higher  value  for 
a.  When  the  boundary  conditions  are  satisfied,  the  temperature 
rise  on  the  back  side  of  the  sample  will  be  greater  than  99%  of 
the  maximum  temperature  rise  after  4  t \ . 

(2)  Relation  between  measured  a  and  Q:  From  the  mathematical 
derivation  in  section  2,  we  know  that  irradiating  the  same  sample 
with  lasers  of  different  Q  should  give  the  same  value  of  a  and  the 
degree  to  which  the  boundary  conditions  are  satisfied  during  the 
test  may  be  judged  from  this.  With  multiple  hole  carbon  as 
a  sample  for  instancy  at  1420°C  we  have  irradiated  with  3  different 
values  of  Q.  The  t \  and  a  obtained  are  all  the  same.  (see  figure 
4).  The  capacitances  of  the  3  experimental  lasers  are  all  20C0  uf. 
The  voltage  values  are  respectively,  1300V,  1400V  and  1500V. 


- —MM  g 


Figure  4.  Experimental 
curve  to  verify  the 
independence  of  a,  from  Q. 

Key:  1 — temperature; 

2 — time. 
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Figure  5.  The  optimum 
thickness  range  of  pure 
iron  at  various  temperatures. 

Key:  1 — sample  thickness 
(mm);  2 — experimental  tempera¬ 
ture  (K). 


2.  Selection  of  sample  optimum  thickness:  To  better  satisfy 
the  boundary  conditions  at  testing  time,  we  must  select  the  optimum 
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sample  thickness.  Too  thick  a  sample  will  increase  the  time  for 
the  temperature  rise  to  reactv-thermal  equilibrium,  increase  the 
heat  loss  and  lead  to  a  larger  value  for  ct.  Too  thin  a  sample 
will  not  satisfy  the  boundary  condition  that  the  laser  pulse  time 
t  should  be  much  less  than  the  time  for  the  sample  to  reach  thermal 
equilibrium,  and  lead  to  a  smaller  value  for  n  .  Now  by  using  pure 
iron  samples  as  examples,  we  shall  calculate  the  range  of  the 
optimum  thickness. 


! 


I 
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(1)  Calculation  indicates  that  when  the  ratio  of  t  to  the 

_2 

characteristic  time  t  Is  less  than  10  ,  the  effect  of  the  time 

c  ’ 

limitation  on  pulse  length  may  be  neglected  [8]  namely 


r/t. 


t 

( L/n)1  •  ®- 


<  io-\ 


(14) 


From  the  above  equation,  for  pure  iron  samples  at  experimental 
temperatures  600  K,  800  K  and  1000  K,  we  calculate  the  optimum 
thickness  at  >1.5,  1.3  and  0.98mm,  respectively. 


(2)  Thermal  diffusivity  calculation  indicates  that  the  radiation 
heat  loss  value  is  determined  by  the  radiation  parameter  y.  When 
y  —  *xi<r*  ,  the  error  introduced  by  radiation  heat  loss  is 
<!#'•'  .  Also 


y  —y.  +  CL/'Yy,* 

(15) 

y»  -  4®  ■  •  L, 

(16) 

« 

y,  —  4«y*r  •  ru  ’  r , 

(17) 

then 

y  1  +■  (  L/r)}.  (18) 

Where  a  is  the  Stef fan- Bolt zmann  constant,  r  is  the  sample 
radius,  s  *  and  cr  are  respectively  sample  axial  and  radial  emissivity 
If  we  are  to  satisfy  y  —  5x10“’  »  then  at  testing  temperatures  of 
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1000  K,  800  K  and  600  K,  the  sample  optimum  thicknesses  are  respectively 
<2.3,3.15  and  7mm. 

Hence,  the  optimum  thickness  range  of  pure  iron  at  various 
temperatures  lies  in  the  shaded  region  in  Figure  5. 


5 .  EXPERIMENTAL  RESULTS 

The  following  are  the  experimental  results  for  3  types  of 
samples  representative  of  the  large  number  of  samples  measured 
with  our  apparatus: 

1.  Armco  pure  iron:  is  the  standard  sample  for  measuring  the  /15 
thermal  diffusivity  of  metals,  and  is  supplied  by  the  Metal  Institute 
of  Academia  Sinica.  The  measured  results  agree  with  the  recommended 
values  of  T.P.R.C.  [9]  and  those  of  Powell  [10]  (see  Table  1  and 
Figure  6 ) . 


Table  1.  Percentage  difference  between  measured  values  of 

ARMCO  pure  iron  and  the  recommended  values  of  T.P.R.C. 
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Figure  6.  Thermal  diffusivity 
curve  for  ARMCO  pure  iron. 

Key:  1 — (calorie/sec  cm  deg) 

2 —  •  our  measured  value 

x  T.P.F..C.  recommended 
value  [9] 

A  Powell  recommended  value 
[10]  j 

3 —  temperature. 


Figure  8.  Thermal  diffusivity 
curve  for  a-AljO^- 

Key:  1 —  (cal/sec*cm  deg); 

2 —  Our  measured  valve; 

1,11  -  the  measured  valves 
from  2  thermal  diffusivity 
apparatus  developed  at  our 
institute;  2 

3 —  Temperature  °C(10  ). 


9  10  n  1;  n  h  is 
j  a*(  *  wk> 


Figure  7.  Thermal  diffusivi 
curve  for  aA^O^* 

Key:  1 — (cal/sec  cm  deg); 

2 —  o  our  measured  value 

x  T.P.R.C.  recommended 
value ; 

3 —  temperature . 


Figure  9.  Thermal  diffusivity 
curve  of  plasma  sprayed  Cr^O^ 
coating. 

Key:  1—  CM2/sec;  2— Sample; 

3 —  Temperature  (X10‘~°C)< 
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(2}  o-AIjO,:  :  is  the  standard  sample  for  thermal  diffusivity 

measurement  of  non-metallic  materials  and  is  prepared  by  our 
Institute.  The  density  is  3.96  gm/cm^,  purity  ^99.5%.  The 
measured  results  agree  with  the  T.P.R.C.  recommended  values  and 
other  thermal  diffusivity  measurements  made  at  our  Institute 
[11,12]  (see  Figure  7,  Figure  8,  and  Table  2). 

Table  2.  Experimental  values  of  *  for  «-ai.o,  and  the 

percentage  difference  from  T.P.R.C.  recommended 


values . 

* 

Temperature  (k)  j  iuoo  nuu  _  1200 

DUO 

1  H/U 

Percentage  1  +i.»  |  +•»•'  |  +5-7 

1 

5.i 

r 

5 .7 

difference  (%) 

(3)  Plasma  sprayed  with  Cr203  coating:  Sample  air  orifice 
rate  and  density  are  respectively  5%  and  5.7  grn/cm^.  The  measured 
values  are  shown  in  Figure  9. 

5.  ERROR  ANALYSIS 

The  weights  of  total  error  in  measurement  of  each  single 
term  errors  are  all  different  for  different  samples  at  different 
temperatures.  We  shall  now  analyze  the  error  in  the  measurement 
for  pure  iron  at  971K. 

1.  The  measurement  error  introduced  by  secondary  apparatus: 

(1)  Error  |«,|:  is  measuring  the  thickness  L:  l  —  2mm 
Thickness  is  measured  with  a  0  level  one  thousandth  meter  bar  with 
maximum  error  ±  in  ,  i.e.,  AL  —  0.002  mm  hence, 

|s, |  -  2 1 A/./L  |  -  ±0.2%. 

(2)  Error  |sil  in  measuring  :  It  includes  5  terms  namely 


(1)  Error  in  the  quartz  crystal  clock  10”^  sec.  It  may 
be  neglected  compared  to  the  value  of  110  ms  for  t  1/2; 

(2)  A  to  D  converter  error  ±0.5%; 

(3)  DC  amplifier  error,  including  the  DC  enhanced  precision, 
enhanced  stability  and  DC  enhanced  linearity  are  respectively  ±u.i%, 
with  a  total  of  ±0.3%; 

(4)  The  computer  stores  one  data  point  every  5  ms,  intro¬ 
ducing  in  t  1/2  an  error  of  2.5  ms.  Since  — 110  ms,  the  error 
introduced  is  ±2.27%; 

(5)  Light  sensitive  resistor  and  differential  amplifier 

time  constants  are  of  the  order  of  micro-seconds.  The  error  intro¬ 
duced  in  t  1/2  may  be  neglected. 

The  maximum  relative  error  in  the  above  5  terms  |6,|  —  ±3.1%. 

2.  Error  introduced  because  the  boundary  conditions  of 
Equation  (12)  cannot  be  completely  satisfied: 

(1)  Error  introduced  by  the  effect  of  time-limited  pulse 

l*»l:  calculation  indicates  that  when  r/r,<l(Tl  ,  this  error 

may  be  neglected.  From  measurement,  r*>0.6ms  ,  .  r,  may  be 

found  from  the  equation  below 

t,  **■  (  L,'*)1  •  a"1  ■“  250ms 

Thus  r  <  10-1  ,  therefore  i«Jj|  0 

t. 

(2)  Error  introduced  by  sample  heat  loss  l««l:  since  the 
experiment  is  carried  out  in  vacuum,  the  heat  loss  due  to  air 
convection  and  conduction  may  be  neglected.  However,  radiation  heat 
loss  Increases  with  the  rise  of  temperature  and  may  be  calculated 
from  Equation  (18),  where  Tt  —  971K.A  'i.o«os  cal/sec  cm 

X2,  <  0.3,  L  —  0.2  cm,  r  —  0.5  cm  .  Hence,  y  —  5.4x10"'.  .  From 

the  relationship  [8]  between  y  and  <*//,  the  constant  term 
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in  Equation  (12)  may  be  found  as  0.137.  Since  the  constant  term 
used  in  Equation  (12)  is  0.133,  the  error  introduced  |e,|  =»  0.73%. 

(3)  Error  introduced  by  the  non-uniform  heating  of  the  laser 
bean  |*,|:  it  is  difficult  to  measure  accurately  the  uniformity 
of  the  beam.  Analyzing  the  photograph  we  took  of  the  beam,  the 
beam  energy  is  quite  uniform.  Besides,  smoothing  of  the  curve 

is  considered  in  the  computer  program.  According  to  the  estimate 

of  Taylor  [2]  ,  the  error  in  this  part  |6j  <  ±i%. 

In  the  5  error  terms  above,  with  the  exception  that  the 
radiation  heat  loss  is  a  positive  error,  the  others  are  all 
positive  and  negative  errors.  Treating  them  by  using  the  root 
mean  square  error,  we  get  the  total  error  as  etota]_  *“  “ 3.2 — 1-3.3%. 

Similarly  sample  at  i443°K  has  a  radiation  error 

l««l  —  5%,  .  With  the  other  terms  respectively  !e,|  ±0.32%, 

Uil  ±2.2%,  Ifiil  —  0,  1 6,1  a*  ±1%.  The  total  root  mean  square  error  is 


ea  —  —2.44%  ~  4-5.6%. 


The  root  mean  square  errors  of  the  2  standard  samples,  and 
the  measured  values  and  the  recommended  values  from  foreign  sources 
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for  the  percentage  errors  are  shown  in  Table  3. 


Table  3  shows  that  compared  to  the  recommended  values  from 
foreign  sources,  the  positive  error  in  the  measured  value  of  \ 
and  its  magnitude  is  in  agreement  with  the  result  of  our  error 
analysis . 


Table  3.  Measurement  error  in  the  2  standard  samples. 


T.P.R.C.  recommended  Percentage  error  between  Our  measured 

error  [9]  our  measured  value  and  error 


T.P.R.C.  recommended 
value 


Pure  iron 

Pure  iron 

Pure  iron 

(97  IK)  j 

a-Al,0,(l448K) 

(971K) 

a-Al.O, 

(97IK.; 

u-Al.O,  (1448K) 

±6— 1U% 

2.6% 

5. 7-5,  | 

-  J.2%— 4-  t. 

-2.44%- 

4-  4 .6% 

7.  RESULTS  AND  DISCUSSION 

1.  After  our  apparatus  is  controlled  by  the  computer,  the 
precision  of  measurement  repeatability  and  speed  of  measurement 
degree  of  automation  have  all  pronounced  improved. 

2.  Since  the  emmissivity  of  the  non-metal  «*  is  large, 

and  its  thermal  diffusivity  1  small,  therefore  at  the  same  measure¬ 
ment  temperature  the  error  introduced  by  heat  loss  is  larger  than 
that  of  the  metallic  material.  There  are  two  ways  to  reduce  this 
error:  deposit  carbonized  tungsten  with  small  8  on  the  sample  or 
correct  by  using  correction  formula  for  radiation  heat  loss. 

3.  Since  the  t \  for  metal  is  smaller,  the  error  introduced 
in  t ?  when  the  computer  takes  in  data  is  greater  than  that  for 
non-metallic  material.  The  way  to  reduce  this  error  is  to  increase 
the  data  sampling  rate  of  the  computer,  or  to  increase  the  sample 
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thickness  without  increasing  significantly  the  radiation  heat 
loss,  i.e.  to  use  the  upper  limit  value  of  the  best  thickness. 
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EXPERIMENTAL  INVESTIGATION  OP  SIMULATING 
IMPINGEMENT  COOLING  OF  CONCAVE  SURFACES 


OF  TURBINE  AIRFOILS 

Cheng-Ji-rui  Wane  Bao-guan 
(Nanjing  Institute  of  Aeronautics  and  Astronautics) 


Abstract 

Presented  in  this  pap«'r  arc  the  experimental  results  of  impingement  cooling  of 
two  uliuninuim  targets  of  different  sizer  with  six  impingement,  tubes.  The  target 
simulating  the  leading  edge  of  a  turbine  airfoil  is  preheated  and  cooled  by  jets  of 
different  impingement  tubes  at  various  flow  rates  and  distances.  The  lumped-capacity 
method  is  used  to  determine  tiie  cooling  rate  of  the  aluminium  target  and  the  mean 
heat  transfer  coefficients  which  are  in  turn  correlated  in  dimensionless  parameters. 
Satisfactory  results  are  obtained  when  Nru  numbers  calculated  from  the  recommended 
expression  are  compared  with  those  by  experiments.  The  recommended  expression  is 
further  compared  with  the  correlated  formula  of  Arizona  State  University  by  Metzgor 
and  those  of  University  of  Cincinnati  by  Ravuri  and  Tabakoff  together  with  the 
experimental  data,  and  it  is  found  that  the  expression  correlated  by  this  paper  is  more 
or  less  reasonable. 


FORWARD 

For  an  aeronautical  jet  engine,  the  thrust-to-weight  ratio 
is  an  important  index  for  judging  the  design  standard.  In  order 
to  improve  the  thrust-to-weight  ratio,  it  is  necessary  to  increase 
the  front  inlet  temperature  of  the  turbire.  Because  the  blade  material 
cannot  endure  too  high  a  temperature,  it  is  often  necessary 
to  employ  cooling  techniques  so  that  the  temperature  in  front 
of  the  turbine  can  be  raised  as  far  as  possible  while  still 
staying  within  a  temperature  endurable  by  the  blade  material. 

Air  convection  cooling  has  been  the  basic  method  for  cooling 
the  aeronautical  combustion  turbine  blades,,  but  for  some  severely 
heated  regions  on  the  blade  (e.g.  the  front  and  rear  edge  of  the 
blade)  convection  cooling  along  is  not  enough.  Other  cooling 
techniques  must  be  incorporated  for  local  cooling.  The  usual 
f 

This  paper  was  presented  at  the  2nd  National  Engineering  Thermal 

Physics  Conference  at  Hangzhon  in  November,  1978. 


method  is  to  use  local  air  film  cooling  (external  cooling)  or 
by  using  air  jets  for  impingement  cooling  (internal  cooling) 
directly  on  the  front  edge  of  the  blade,  or  by  using  both. 

Using  impingement  air  jet  to  reinforce  the  local  cooling  of 
the  leading  edge  is  a  structurally  simple,  practical  and  easy 
method . 

The  theory  of  jets  and  experimental  study  has  long  been 
discussed,  but  the  discussions  are  usually  limited  to  free  jet 
or  wall  jet  and  not  about  impingement  jets  normal  to  flat  plates 
or  concave  surfaces.  In  some  papers,  only  the  flow  characteristics 
of  flat  plate  impingement  jets  are  supplied  while  the  heat  exchange 
problems  of  impingement  jets  on  plane  or  concave  surfaces  are 
not  discussed. 

From  the  end  of  the  sixties  to  the  beginning  of  the  seventies, 
in  order  to  apply  impingement  jets  to  the  cooling  of  the  leading 
edge  of  turbine  airfoil  some  research  was  carried  out  [5-12]  on 
the  impingement  cooling  of  concave  surfaces,  the  majority  of  which 
was  basically  experimental  studies.  Among  these  experimental 
studies,  the  research  works  of  D.E.  Metzger  of  the  University 
of  Arizona  and  those  of  W.  Tabakoff  of  the  University  of  Cincin¬ 
nati  in  the  U.S.  are  more  outstanding. 

Metzger  did  some  experimental  studies  on  the  impingement 
cooling  of  a  row  of  circular  jets  on  concave  surfaces  [6,7].  He 
analyzed  the  effect  of  the  parameters  Re,  'Jd„  l/d;  and 
on  the  Stanton  number  St  (for  the  meaning  of  symbols,  see  the 
symbol  explanation)  and  finally  proposed  the  empirical  formula  [6] 

ScmaiRcul7(Ww-0.35S  (1) 

The  range  of  applicability  is  4.65  <  lib  <  iv6,  1.67  <«•„/</„<  6.67,  i  isn  <  Kc  <  6>oo. 
He  did  not  include  in  eq.  (1)  quantitatively  the  effect  of  cjd„ 
and  z „/i  ,  but  only  used  1/b  to  implement  indirectly  the 

effect  of  cn/dn  and  to  eliminate  the  effect  of  Zn/b  by  taking  the 


/16 


maximum  value  of  St  number. 


Tabakoff,  et  al. ,  also  carried  out  a  more  thorough  analysis 
and  study  [9-12]  over  a  period  of  time  on  the  impingement  cooling 
of  semi-spherical  concave  surfaces.  In  Reference  [10]  he  proposed 

Nu,„  -  cRei-fV'X  cjd,  HF if)  ,  ?  v 


as  the  standard  relationship  for  impingement  cooling  of  semi- 
spherical  concave  surfaces  by  multiple  array  of  circular  jets. 
In  Reference  [11]  he  proposed 

Nu,  =  O.'j  1 1  3Rep  l5(y,,'D)0  '1j( c,,  d,)~  '  i:*(  1  5 


as  the  standard  relationship  for  the  impingement  cooling  of 
semi-spherical  concave  surface  by  single  array  of  circular  jets. 
In  the  above  2  equations, 


=  Vpd./fi 


(M 


RcP  “  VpW'u 


(5) 


where  W  is  the  length  of  the  curved  surface  of  the  impinged 
concave  plate,  i.e. 


iv 


*Dj  2 


(6) 


In  eq.  (2 ) , 

of  Z  /d  : 
n  n 


the  coefficient  c  and  the  index  q  are  both  functions 


C  -  l.36txp(-VzjJJ  -  0.002(Z,/d ,)  -  0.153.  (ZJd,)  '7) 

\.\^{Zjd,)  -  3.23  (8) 

In  the  formulae  proposed  in  Reference  [10,11],  the  diameter 
of  the  jet  orifice  or  the  length  of  the  curved  surface  W  of  the 


impinged  concave  surface  are  taken  respectively  as  the  formative 

dimension  of  Re  and  Mu,  while  the  parameters  such  as  the  number 

of  jet  orifices  n  and  the  center  distance  c  have  not  been  taken 

n 

into  consideration.  This  seems  to  be  a  lack  of  completeness. 

In  this  paper,  we  synthesized  the  results  of  the  papers  men¬ 
tioned  above  and  simulated  the  impingement  cooling  of  the  leading 
edge  of  the  airfoil  by  cooling  air  by  impinging  against  the 
semi-spherical  concave  surface  with  a  single  array  of  circular 
air  jets.  The  experimental  set  up  is  basically  based  on  the 
University  of  Arizona  method  while  the  data  are  treated  orimarily 
according  to  the  University  of  Cincinnati  method. 


EXPERIMENT  AND  MEASUREMENT  SET  UP  AND  THE  EXPERIMENTAL  PROCEDURE 

Figure  1  shows  the  diagram  of  the  air  passage  system  in  the 
experimental  set  up.  The  air  stream  from  the  pressure  stilling 
chamber  passes  through  the  regulator  valve  4,  the  rotor  flow 
rate  meter  7,  the  air  storage  chamber  10  and  impingement  tube 
17  and  impinges  on  the  test  sample  16  to  simulate  cooling  of  the 
leading  edge  of  the  airfoil.  Finally  the  cooling  air  stream 
is  expelled  to  the  surrounding  environment.  The  rotor  flow  rate 
meter  and  the  nickel-chromium-copper  thermal  couple  have  all 
been  calibrated. 


Figure  1.  Diagram  of  the  experimental  set  up. 

Key:l — source  of  air;  2 — pressure-stilling  chamber;  3 — pressure 
gauge;  4 — regulating  valve;  5 — heating  oven;  6 — pressure  gauge; 
7__r*otor  flow  rate  meter;  8 — digital  voltmeter;  9 — soft  turbine 
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Key  to  Fig.  1  cont.:  10 — air  storage  chamber;  11 — thermometer; 
12 — pressure  gauge;  13 — work  bench;  14 — thermal  couple; 

15 — insulated  chamber;  16 — test  sample;  17 — impinge  tube; 
18 — coordinate  regulating  stand. 


Figure  2  shows  the  diagram  of  the  test  piece,  the  impingement 
tube  and  their  set  up.  The  semi-circular  concave  surface  on 
the  test  piece  is  impinged  directly  by  the  cooling  air  stream. 

Our  experiment  used  2  test  pieces  with  semi-circular  concave 
surfaces  of  different  diameters  D:  D=15mm  and  D=llmm.  The 
test  piece  is  heated  first  in  a  small  electric  heater.  Its  tempera¬ 
ture  is  determined  by  the  nickel-chromium-copper  thermal  couple 
buried  in  the  test  piece  with  a  PZ8  digital  DC  voltmeter.  After 
being  heated,  the  test  piece  is  removed  from  the  electric 
heater  and  placed  into  the  cavity  of  an  insulated  box  made  from 
Tung  wood.  The  cover  of  the  insulated  box  is  then  closed  with 
only  the  semi-circular  concave  surface  of  the  test  piece  exposed 
to  the  jet  of  the  impingement  tube  and  cooled  by  it.  The  tempera¬ 
ture  change  of  the  test  piece  is  displayed  on  the  PZ8  digital 
voltmeter  through  the  thermal  couple.  After  the  impingement  air 
stream  flow  rate  is  stabilized  and  fixed  at  a  required  volume 
by  adjusting  the  regulating  valve,  the  Instantaneous  test  piece 
temperatures  are  recorded  with  the  timer  and  the  digital  voltmeter 
so  that  the  variational  relationship  between  the  instantaneous 
temperature  of  the  piece  and  time  t  is  obtained,  from  which  the 
heat  exchange  coefficient  o .of  the  impingement  jet  on  the  test 
piece  averaged  over  the  semi-circular  concave  surface  is  found. 

This  a  value  corresponds  to  the  average  heat  exchange  coefficient 
at  a  certain  flow  rate.  Thus  different  a  values  for  different  flow 
rates  are  obtained. 

The  width  d  of  the  impingement  tube  is  1/2  of  the  diameter 
of  the  semi-circular  concave  surface  of  the  corresponding  test 
piece,  i.e.  d/D-1/2.  Hence  there  are  2  different  widths  for  the 
impingement  tube:  d*7.5  and  d*5.5mm.  The  orifice  diameter  of 
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the  impingement  tube  is  determined  by  a  tool  microscope  and 
an  average  value  d^  is  then  obtained.  For  an  impingement  tube 
of  the  same  width,  there  are  3  different  center  distances  Cn 
i.e.,  cjdm  =  1.67, 3.33  and  6.67.  The  impingement  tube  is  fixed  on 
a  tri-directional  coordinate  regulator  stand  18  so  that  the  dis¬ 
tance  from  the  impingement  tube  to  the  semi-circular  concave 
surface  of  the  test  piece  may  be  adjusted  at  will.  Each  impinge¬ 
ment  tube  is  first  adjusted  to  be  at  the  given  distance  Z^. 

Then  a  sec  of  average  heat  exchange  coefficients  corresponding 
to  a  given  flow  rate  is  measured  and  their  average  value  is 
taken  to  be  the  average  heat  exchange  coefficient  at  that  flow 
rate  to  reduce  accidental  error.  This  is  done  for  different 
flow  rates.  Then  the  above  experiment  is  repeated  for  another 
Zn  value.  In  our  experiment  there  are  2  test  pieces  and  for 
each  test  piece,  there  are  3  impingement  tubes.  Five  Zn 
values  are  used  for  each  tube  and  for  each  Z  value,  exDeriments 
are  carried  out  for  five  different  flow  rates.  Hence,  one 
empirical  formula  may  be  found  for  each  Zn  value,  for  a  total 
of  30  empirical  formulae,  (for  detail,  see  Table  1). 


TREATMENT  OF  EXPERIMENTAL  DATA 

The  test  pieces  are  made  from  pure  aluminum.  As  the  thermal 

conductivity  of  aluminum  blocks  is  large  and  the  test  piece 

dimension  is  small,  the  average  heat  exchange  coefficient  obtained  /168 

2 

from  this  experiment  is  mainly  between  170  to  1200  Kcal/m  •  time  . 
degree.  Preliminary  estimation  gives  Biot  number  8.  —  0.025  ~  0.LZ5  > 
averaging  about  0.1.  The  non-uniformity  in  temperature  is  no  more 
than  2.5%.  The  characteristic  dimension  of  the  aluminum  block 
is  obtained  from  the  ratio  of  the  aluminum  block  volume  to  the 
cooled  heat  exchange  surface.  We  may  assume  that  the  spatial 
temperature  distribution  of  the  aluminum  block  is  basically 
uniform  so  that  the  cooling  rate  of  the  aluminum  block  may  be 
calculated  from  the  concentrated  total  thermal  capacity  method. 

The  thermal  equilibrium  equation  of  the  test  piece  during  cooling  is 

—  G.n.c ,l  •  dT /dt  —  «F(  T  —  r„)  /  0  \ 


Integrating  the  above  equation,  we  get 


In [(7,  -  T-)/(r  -  T »)]  -  <*F(t  -  tdiC'ALCAL 


(10) 


where  T^  is  the  initial  temperature  at  the  initial  instantaneous 
time  t.^,  Gal  and  cAL  are  respectively  the  weight  and  specific 
heat  of  the  aluminum  block.  Let 


c  =  G.tt  c  n.i  F  » 

and  take  into  account  of  heat  loss  by  leakage,  the  above  equation 
may  be  re-written  as 


[c/O  ~  O ] In [ ( x,  -  T.)/(T  -  r«)]  - 


(11) 


The  leakage  heat  loss  is  measurable:  Place  the  aluminum  block 
into  the  insulated  box  after  heating,  block  the  semi-circular 
concave  surface  with  Tung  wood  and  let  it  cool  by  itself.  The 
required  time  interval  for  T.^  to  decrease  to  T  is  measured  and 
the  leakage  heat  loss  is  estimated  from  it.  Experiment  indicates 
that  the  leakage  heat  loss  is  very  small,  about  3%  of  the  heat 
exchange  coefficient  at  the  minimum  flow  rate.  It  may  be  neglected 
within  the  accuracy  of  our  experiment,  when  the  impingement  jet 
flow  rate  has  stabilized,  time  measurement  with  the  second  timer 
is  started  (take  t^=0)  and  simultaneously  the  instantaneous  tempera¬ 
ture  Ti  of  the  aluminum  block  is  also  recorded.  Then  the  aluminum 
temperatures  at  a  sequence  of  pre-arranged  time  intervals  are 
recorded  and  a  set  of  heat  exchange  coefficients  for  a  given 
value  of  impingement  jet  flow  rate  is  calculated  according  to 
Eq.  (11)  from  which  the  arithmetic  average  value  is  taken  to  be 
the  average  heat  exchange  coefficient  at  this  flow  rate.  The  Nu 
number  is  then  calculated 


Nu  =*■  “  a  •  2b j  A, 


(12) 


where  d  is  the  equivalent  diameter  when  the  single  array  circular 

a 
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b 


hole  Is  converted  Into  the  equivalent  2-dimensional  jet  channel 
and  b  is  the  width  of  the  equivalent  2-dimensional  jet  channel, 

i.  e. 


b  “  n(  iti  4  )d l  „ 


(13) 


The  Re  number  of  the  jet  is 


KC  =  V j,iir  n  -*m{G  •  lb)  /.'<*>' 


(14) 


The  volume  flow  rate  measured  by  the  rotor  flow  rate  meter 

p 

is  standardized  by  the  standard  state  PQ=10332  kg/m  ,  Tq=293K. 

This  volume  flow  rate  value  should  be  corrected  according  to  the 
air  stream  state  during  experiment  and  converted  to  the  weight 
flow  rate  G  according  to  the  current  air  stream  pressure  and 
temperature  with  the  following  equation: 


G  —  0.2029 g„  v  P  7 


(15) 


where  P  is  the  absolute  pressure  (kg/nr  )  of  the  air  storage  chamber, 
T  the  absolute  temperature  (K)  of  the  air  storage  chamber  and 
Qq  the  volume  flow  rate  at  STP. 

According  to  similarity  theory,  the  heat  exchange  coefficient 
may  be  expressed  with  the  following  equation 


Nu  —  /(Ke,  Pr,  if»),  <pi,  T  »  '  * ' ) 


(16) 


For  air,  the  Prandtl  number  Pr  basically  does  not  vary  very  much. 
We  may  take  its  value  as  fixed.  <p\,  >pt,  <jv  •  •  are  geometrical 
similarity  standards.  According  to  experimental  analysis,  we  take 
the  following  similarity  standards: 

lf),  —  ZJb  j  the  relative  distance  between  the  test  piece  and 
jet  orifice; 


36 


>f>,  “i,.  dn  ,  the  relative  center  distance  of  the  jet  orifices; 

:f-  =  F ■  )  ,  ratio  of  the  semi-circular  concave  surface  area 

to  the  total  cross-sectional  area  of  the  jet  orifice,  which 
reflects  the  effect  of  the  diameter  of  the  concave  surface  on 
a.  Eq.  (16)  may  be  written  as 

Nu  =  .l,Re '*(<  „/</,;  '•(  7.  ,  b)  '.(  F  t)A' 

(17) 

For  a  given  impingement  tube  and  test  piece,  cn/cin  and  ’s 

values  are  all  fixed.  For  a  given  Z^/b  Eq.  (17)  may  be  re-written 
as 

Nu  =  ./0Ke-*‘ 

(18)  71-69 


where 


=  Ax(c,„J,)'<ZJb)HF  7)-'« 


(19) 


EXPERIMENTAL  RESULTS  AND  ANALYSIS 

Table  1  shows  the  empirical  formulae  obtained  from  the  data 
of  6  impingement  tubes  at  5  different  Z^/b  values  for  5  different 
flow  rates  in  the  experiment,  after  re-ordered  into  the  form  of 
Nu=ARem.  Because  of  the  effects  of  different  experimental  condi¬ 
tions  and  measurement  errors,  etc.,  the  A  and  m  values  of  each 
formula  will  not  be  the  same. 

Figures  3  through  7  show  the  comparison  between  experimental 
data  and  the  empirical  formula  obtained  from  them  for  the  6  impinge¬ 
ment  tubes  at  a  fixed  Zn/b  value  (3  types  with  cJdK  —  i.67,  3U3,  6.67  ) 

Details  about  the  meaning  of  the  symbols  in  the  figures  may  be 
found  in  Table  1.  Generally  speaking,  Nu  number  decreases  ».  j.th 
the  increase  of  zn/b  >  but  Zn/b  ’ s  value  varies  from  3  to  80 
while  cn/dn  from  1.67  to  6.67.  Hence,  relatively  speaking,  the 

effect  of  Z  /b  on  the  Nu  number  is  not  as-  large  as  that  of  c  /d  . 

n  n  n 
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Table  1.  Detailed  table  of  Figures  3-6 


In  Figure  7,  with  the  exception  of  line  16  labeled  with  the 
experimental  value  of  Zn/b=3,  the  other  5  lines  all  are  the 
experimental  data  at  Zn/b=80. 

From  Figures  3  through  7,  we  can  observe  the  effect  of 
c  /d„  on  Nu:  the  Mu  number  decreases  significantly  with  the 
increase  of  cn/dn.  Under  the  condition  of  equal  cn/dn  value, 
the  Nu  numbers  of  test  pieces  are  relatively  close.  The  Nu 
number  of  the  test  piece  with  D-llmm  is  slightly  higher  than  that 
of  the  test  piece  with  D-15mrrt.  From  this  it  can  be  seen  that 
Nu  number  is  related  to  the  heat  exchange  surface  area  of  the  test 
piece,  and  hence,  related  to  the  radius  of  the  semi-circular 
concave  surface  of  the  simulating  airfoil’s  leading  edge. 

Metzger  [6]  has  pointed  out  that  the  average  heat  exchange 
coefficient  of  the  impingement  jet  varies  with  Z^/b;  for  a  fixed 
value  of  cn/dn,  a  maximum  heat  exchange  coefficient  is  found 
for  some  Zn/b  value.  After  passing  through  the  value,  a 
decreases  with  increment  of  Z  /b. 


R* 
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The  value  of  Zn/b  producing  the  maximum  heat  exchange  coeffi-  /171 
cient  is  related  to  cn/dn*  When  cJd .  —  1-67  ,  the  maximum 

value  of  heat  exchange  coefficient  occurs  at  l  <  ZJb  <3  ;  when 

cjdn  —  3.33  ,  it  occurs  at  2<Z„/£<5  and  when  —  6.67  » 

it  occurs  at  5s$z„/*<n  .  in  our  experiment,  only  for  the 

test  piece  with  D=15mm  did  it  reach  Z^/b-3  once  when  cjdn  —  1.67 
With  the  exception  of  this,  it  is  very  difficult  to  achieve  Zn/b*3 
in  the  experiment;  therefore  it  is  not  possible  to  determine  the 
Zn/b  value  through  experiment  that  produces  the  maximum  heat 
exchange  coefficient.  The  above  is  the  empirical  formulae  obtained 
for  each  impingement  tube  at  a  given  value  of  Z^/b  in  the  experiment: 

Nu  —  ARcn, 

The  values  of  A  and  m  vary  for  various  values  of  c  /d„  and  Z  /b- 
hence,  30  empirical  formulae  are  obtained.  They  are  not  convenient 
to  use  although  it  is  easier  to  compare  them  to  experimental  data. 

Hence,  we  took  all  the  data  (150  experimental  points)  obtained 
from  experiment  and  put  them  in  a  matrix  form  according  to  the 
least  square  method  based  on  Eq.  (17)  and  finally  find  the  solution 
through  the  electronic  computer  with  a  program  (which  is  here 
omitted).  The  coefficients  and  indices  in  Eq.  (17)  are  obtained 
respectively  as 

A\  —  0.269,  At  —  0.666,  ,*,--0.401,  ,*,--0.065,  /*,  —  -0.201. 

(20) 


i.  e. 


(21)  /vn 


The  applicable  ranges  of  the  above  equation  are: 

1950  <  Re  <  8900,  1.67  <  ejdn  <  6.67,  3  <  ZJb  <  80,  36  <  Fjj_  <  135.  . 

The  qualitative  temperature  is  taken  as  the  total  temperature 
(air  temperature  of  the  air  storage  chamber)  of  the  impingement 
air  stream  in  front  of  the  orifice.  Prom  the  index  of  Eq.  (21) 
one  can  see  intuitively  that  besides  Re,  the  factor  that  affects 
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the  heat  exchange  coefficient  the  most  is  c  /d„  while  the  effect 
due  to  Zp/b  *s  fcke  least.  To  test  the  agreement  of  the  composite 
relation  (21)  with  experimental  data,  we  compare  the  Nucal  calcul¬ 
ated  from  Eq.  (21)  to  the  corresponding  Nu  measured  from  experi- 

w  Ap 

ment.  The  experimental  data  points  are  shown  with  the  symbols 
in  Table  1,  as  shown  in  Figure  8.  From  Figure  8,  it  can  be  seen 
that  the  maximum  deviation  between  NucgQ  and  Nuexp  is  within 
±17$.  90$  of  the  experimental  points  fall  within  ±8$.  This 

result  is  satisfactory.  The  index  of  Re  in  the  composite  formula 
(17)  is  precisely  the  middle  value  of  the  30m  values  in  the 
empirical  formula  listed  in  Table  1.  The  size  of  the  index  shows 
that  the  flow  state  of  the  jet  stream  is  close  to  turbulent  flow. 

If  the  Eq.  (21)  is  to  be  extended  to  the  impingement  cooling 
by  other  fluid,  it  can  be  done  by  the  traditional  method  by  adding 

1  /  "3 

on  the  LHS  of  the  equation  a  factor  Pr  ' J  as  a  correction.  The 
coefficient  A^  should  be  changed  to  0.326  namely 

Nu  -  0.326Re4WPr (21a) 

COMPARISON  EQ. (21)  OF  THIS  PAPER  TO  REFERENCE  [6]  AND  [11] 


Metzger  pointed  in  Reference  [6]  that  the  average  heat  exchange 
coefficient  varies  with  Zn/b  .  Yet  when  the  data  are  organized, 
only  the  largest  heat  exchange  coefficient  value  for  each  Zn/b 
value  and  the  effect  of  Zn/b  is  not  incorporated  with  the 
formula.  The  formula  realizes  the  effects  of  cn/dn  and  F/f 
indirectly  through  the  parameter  1/b.  In  practical  applications, 
the  value  of  Zn/b  that  gives  the  maximum  heat  exchange  coefficient 
is  usually  unknown  and  the  effects  of  cn/dn  and  F/f  cannot  be 
expressed  in  the  formula. 

In  order  to  make  it  easier  to  compare  Eq.  (1)  to  the  Eq.  (21) 
recommended  In  our  paper  and  to  the  experimental  data,  it  is 
necessary  to  transform  the  equation  as  follows: 


Nu...R*’ *(//*)* M  -  0.355.  RePr 


Nu„.„Re,,7(//A)»»  —  0.355  RePr 


(22) 


Combining  the  case  of  D^lSmm  in  our  experiment,  in  air  Pr»0.7 
in  Eq.  (22) 

/  -  (1/2)(1/2)(*D)  -  1 1.76mm, 

b  is  related  to  the  impingement  tube  used.  Substituting  the  b 
values  corresponding  to  the  3  cn/dn  values  into  Eq.  (22),  the 
Metzger  equation  may  be  changed  to: 


(a)  ^  cjd,  —  1.67  ft,  Num„  —  0.0577Re°-71 

(b)  ^  cjdn  —  3.33  ft,  Num„  —  0.04025Rea-7’ 

(c)  ^  cjd,  —  6.67  ft,  Numai  —  0.0307Reu-73 


(23) 


Figure  9  shows  the  comparison  of  Eq.  (21)  in  this  paper  to 
the  Metzger  Eq.(23).  In  the  diagram,  it  can  be  seen  that: 

When  1 67  t  gq#  (23)  is  close  to  Eq.  (21)  and  to  the 

experimental  data  in  this  paper.  When  — 3.33  ,  Eq.  (23) 

is  slightly  higher  than  our  Eq.  (21)  or  (18)  at  Zn/b»5,  and  when 
cm/da  —  6.67  ,  Eq.  (23)  deviates  more  from  (21)  or  (18). 

This  is  because  the  heat  exchange  coefficient  in  Eq.  (23)  is 
meant  to  be  the  maximum  value.  When  ejd„  — •  1.67  t  the  maximum 

value  occurs  when  Zn/b  is  between  1  and  3 »  hence,  Eq.  (23)  is 
closer  to  Eq.  (21)  at  Zn/b«3.  When  c./rf.  —  6.67  ,  the  maximum 

value  situation  did  not  occur  and  Eq.  (23)  is  higher  than  Eq.  (21) 
Through  comparison,  it  is  recognized  that  it  is  impossible  to 
include  all  the  different  geometrical  conditions  with  a  single 
geometric  parameter  1/b  . 
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Figure  9.  Comparison  of 
Eq.  (21)  and  Eq.  (23)  of 
Reference  [6]. 

Key:  1 — our  Eq.  (21); 

2— Eq.  (23),  Ref.  [6]. 
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Figure  10.  Comparison  of 
Eq.  (21a)  and  Eq.  (3b)  of 
Ref.  [11]. 

Key:  1 — our  Eq.  (21a); 

2— Eq.  (36,  Ref.  [11]. 


The  empirical  Eq.  (3)  proposed  in  Reference  [11]  covers  the 
experimental  range 

cjdn  -  2.5 — 5.0  Z./d .  -  2-6 

For  the  convenience  of  making  comparisons,  we  choose  our  experimental 
data  closer  to  the  experimental  range  of  Eq.  (3): 

D  —  \5mm.ejd.  —  3.33,  Zjh  —  5,  10,  20  and  40  (corresponding  to 


the  cases  of  zjd.  —  1.25, 2.5, 5.0  and  10.0  In  Reference  [11], 

Eq.  (3)  must  also  be  changed  to  the  dimension  2b  adopted  in  our 
paper: 

Nu  -  0M\lMnKW/2b)^*KdJDy  ™(b/dn)'«  x^Zm/by*-'3\cJd.y*M 

(3a) 


Since 


fp  —  nD/2  —  3.14  X  15/2  —  23.55  mm 
J„  —  1.48mm,  b  *=  0.37mm,  F/f  —  68.3 


Substituting  the  above  values  in  Eq.  (3a)  and  simplifying,  we  get 

Nu  -  0.0589R t*KcJ4.Y'-aKZjb)’*x*  (3b) 

Substituting  again  F/f-68.3  into  Eq.  (21),  we  get  after  simplifi¬ 
cation 

•  Nu-0.1153Re»‘«(r./rf.r^(Z./*)-«  (2la) 

Figure  10  shows  the  comparison  between  our  (21a)  and  3b)  of 
Reference  [11]  with  the  experimental  data  in  our  paper  also 
indicated  in  the  diagram.  From  the  diagram  it  can  be  seen  that 
Eq.  (3)  of  Reference  [11]  is  lower  than  Eq.  (21)  recommended  in 
our  paper  and  the  experimental  data,  and  is  close  only  for  the 
cases  of  Zn/b*5  and  10.  This  is  apparently  due  to  the  fact  that 
Eq.  (3)  in  Reference  [11]  is  only  applicable  for  Zji,  —  2—6 
which  corresponds  to  the  case  of  zjb  —  8—24  .  Therefore  when 

Zjb  >  20  ,  Eq.  (3)  is  not  really  applicable. 

Since  the  impingement  jets  introduced  in  Reference  [10]  have 
been  multiple  array  orifices,  it  is  too  different  from  our  experi¬ 
mental  environment  and  no  comparison  can  be  made. 


CONCLUSIONS 


A7i 

Our  major  conclusions  are  as  follows: 

1.  The  major  parameters  affecting  the  impingement  heat 
exchange  coefficient  are  Re,  cjdm,  Zjb  and  F/f.  The  geometric 
parameters  having  the  largest  effect  on  Nu  number  is  cn/dn,  next 
is  Zn/b  and  the  least  is  Zn/b. 

2.  When  the  experimental  data  are  treated  with  the  least 
square  method  through  computer  calculations,  the  standard  relational 
equation  obtained  is 

Nu  -  0.269Re4  “(r./^)-#  *,(Z./i)-,  “,(i?//)-“J5r 

3.  When  the  results  calculated  with  the  above  equation  are 
compared  to  the  experimental  data  in  this  paper,  basic  agreement 
is  obtained.  When  the  above  equation  is  compared  to  those  recom¬ 
mended  in  Reference  [6]  and  [11],  our  equation  is  considered  to 
be  closer  to  the  experimental  data. 

4.  The  equation  proposed  in  our  paper  basically  includes 
a  number  of  factors  affecting  the  heat  exchange  coefficient  of 
impingement  Jet  in  a  composite,  standardized  relational  formula. 

It  has  a  definite  practical  value. 

SYMBOL  EXPLANATION 

At  a  i,  •••  constant 

b  equivalent  2-dimensional  jet 

c  coefficient  <•  -  galcal/f 

c  center  distance  of  impingement 
n  tube  orifice 

CAL  sPecific  heat  of  aluminum 

D  diameter  of  semi-circular 
concave  surface  of  test  piece 


P  absolute  pressure 
Pr  Prandtl  number  iv-p gCt/x 
Q  Volume  flow  rate  at  STP 

(?,  -  10332,  T  -  293K) 

q  index,  see  Eq.  (2)  and  (8) 

Re  Reynold  number  R«  -  »'/*#/**  -  c  •  :*///*« 
Rep  Rep  -  VpW/y 

T  instantaneous  time 


d  impingement  tube  thickness 

d  equivalent  diameter 
e 

d„  jet  orifice  diameter 
n 

F  heat  exchange  surface  area 

f  jet  orifice  total  cross- 
sectional  area 

G  Jet  flow  rate  by  weight 

G -T  weight  of  aluminum  block 
L  (test  piece) 

g  gravitational  acceleration 

L  length  of  test  piece 

1  length  of  semi-circular  arc 
of  the  heat  exchange  surface 
of  test  piece  /  ■»  (i/2)(*c/2) 

1  outer  distance  between  the 
n  2  ends  of  the  impingement 
tube 


m  power  index 
Nu  Nusselt  number 


Nu  <■  a  «  2*/X 


St  Stanton  number  st  -  a/vre,  (36oo) 

SIm,  —  a^/VrefCiSOO) 

T,  initial  temperature  of 
1  test  piece 

T 

»  total  temperature  or  ambient 
temperature  of  jet  in  the 
air  storage  chamber 

t  time 

V  velocity 

W  arc  length  of  semi-circular 
concave  surface  of  the  test 
piece  H/-ao/2 

Z  distance  from  jet  orifice 
n  outlet  to  test  piece 

a  average  heat  exchange  coefficient 

X  thermal  conductivity 

p  density 

y  dynamic  viscosity 


n  number  of  Jet  orifices  on  the 
impingement  tube 
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Abstract 

A  lot  of  theoretical  and  experimental  investigations  in  relation  to  the  subject  have 
been  done,  and  the  results  demonstrate  that  the  present  solution  for  the  one  dimensional 
simultaneous  equations  is  reasonable.  Both  the  widely  current  flow  drug  calculating 
method  and  the  hole  blank  off  experimental  method  which  was  developed  on  the  base 
of  the  above  said  method  have  some  defects. 

This  paper  provides  a  simple  calculating  method  which  utilizes  the  meau  flow 
coefficient.  The  method  is  simple  in  calculation,  clear  in  concept  and  lias  a  wide  range 
>>f  applications.  The  calculating  result  is  in  good  agreement  with  the  solution  for  the 
Hlie  dimensional  simultaneous  equations. 


1.  EXPERIMENTAL  SET  UP,  TEST  PIECE  AND  MEASUREMENT  SCHEME 

This  experiment  was  carried  out  on  a  water  flow-simulation 
experimental  unit.  The  test  piece  is  a  2-dimensional  test  section 
made  of  organic  glass  with  the  2  passages  separated  from  the  flame 
tube  with  a  separator. 
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The  flow  field  of  the  2  passages  is  measured  by  the  total 
static  pressure  difference.  On  each  cross-sectional  area,  total 
pressure  is  measured  at  9  points.  The  flow  rate  coefficient  over 
the  cross-sectional  area  is  standardized  a  priori  with  Venturi 
tube  flow  rate  ;  then  the  flow  rate  through  each  cross-section 
may  be  measured  during  the  actual  experiment.  The  flow  rate 
through  a  small  orifice  is  determined  by  the  difference  between 
the  flow  rates  of  the  passage  cross-section  before  and  after 
the  orifice. 

*This  paper  was  presented  at  the  2nd  National  Engineering  Thermophysics 
Conference  at  Hangzhon  in  November,  1978. 


The  orifice  jet  velocity  and  orifice  flow  rate  coefficient 
are  treated  according  to  the  cross-section  area  parameter  before 
the  orifice. 

2.  EXPERIMENTAL  RESULTS  AND  ANALYSIS 

Figure  1  shows  the  relation  between  the  orifice  flow  rate 
coefficient  c  and  flow  rate  parameter  F  .  This  result  is  very 
close  to  the  data  given  in  most  papers  (e.g.  Reference  [1])  with 
the  c  value  at  large  Fp  slightly  higher.  (Note:  the  internal 
flow  of  the  flaming  tube  in  our  experiment  is  not  zero.) 

Figure  2  shows  the  relation  between  the  orifice  flow  rate 
coefficient  c°  and  area  parameter.  The  curve  given  in  Reference 
[2]  is  also  drawn  in  this  diagram  (Miller  curve).  From  the 
diagram,  it  can  be  seen  that  for  small  Gp  values,  the  experimental 
values  of  the  orifice  flow  rate  coefficient  are  far  lower  than 
the  Miller  values . 

The  variational  relations  between  the  flow  rate  distri¬ 
bution  W^  of  various  test  pieces  and  the  relative  area  of  the 
orifice  of  the  flaming  tube  J.  are  shown  in  Figure  3  and  Figure 
We  only  emphasize  one  point  here:  namely  that  for  a  fixed 
a^^  value,  the  bigger  the  curve  Xmax  »  the  more  curved  is  the 

curve  w,  —  1(a)  .  The  smaller  the  value  of  X _ _ ,  the  more 

w i  «■»  f(A)  approaches  a  straight  line  of  45°. 

Figures  5  and  Figure  6  respectively,  show  the  value  of  the 
combustion  chamber  friction  coefficient  .  It  can  be  seen 

that  in  general,  the  total  pressure  loss  of  the  combustion  chamber 
is  usually  larger  than  the  velocity  head  at  the  inlet  of  the  2 
passages. 

Now  we  shall  discuss  an  experimental  method  introduced  in 
some  references  to  measure  the  flow  rate  distribution  in  a 


combustion  chamber  by  a  orifice-blocking  method  which  is  the 
so-called  orifice  blocking  method.  The  method  is  as  follows: 
the  total  friction  coefficient  5  in  the  combustion  chamber  is 
measured  first.  Then  the  ith  row  of  orifices  is  blocked  and 
the  friction  coefficient  £  is  again  measured.  From  the  equation 
in  Reference  [2] ,  we  should  mention  that  this  method  is 
quite  questionable.  In  fact,  we  can  easily  see  that 
this  method  basically  is  incapable  of  showing  the  difference 
of  the  flow  rate  coefficient  C  corresponding  to  the  different 
orifice  arrays  at  positions  in  front  of  or  behind  the  flaming 
tube.  Experimental  results  verified  this  completely.  From 
Table  1  we  see  that  the  flame  tube  flow  rate  distribution  measured 
by  the  orifice  blockage  method  is  basically  the  same  for  each 
array  of  orifices.  In  other  words,  the  flow  rate  distribution 
measured  by  the  orifice  blockage  method  is  in  fact  distributed 
according  to  area.  It  deviates  quite  severely  from  accurately 
measured  experimental  results  (respectively,  0.013  and  0.504). 

Of  course,  when  the  total  orifice  area  of  the  flaming  tube  is  very 
small,  the  law  of  actual  flow  rate  distribution  is  very  close 
to  be  law  of  distribution  by  area.  Here  the  orifice  block  method 
naturally  agrees  with  the  experimental  value,  but  what  meaning 
does  this  agreement  have? 


Figure  1.  Relation  between 
orifice  flow  rate  coefficient 
and  flow  parameter. 
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Figure  5.  Relation  between 
friction  coefficient  £  and 
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Figure  6.  Comparison  between 
calculated  values  and  experimental 
values  of  eal  «majt-6.523). 
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Table  1.  Experimental  result  of  orifice  blockage  method 

0»,  =  0.6i,  *„,„=•  6.523) 


\  parameter 

sxper^v 
nental  'v 
aethod 

€ 

r 

AW1,2 

AW11,12 

to  orifice  is 
slocked 

1  1.13 

0.0420 

0.504 

Lst  and  2nd 
row  of  orifices 
slocked 

— 

1.18 

0.0214 

— 

Llth  &  12th  row 
sf  orifices 
blocked 

— 

1.16 

— 

0.013 

METHOD  OF  CALCULATING  THE  COMBUSTION  CHAMBER  FLOW  RATE  DISTRI¬ 
BUTION  AND  THE  AVERAGE  FLOW  RATE  COEFFICIENT  OF  TOTAL  PRESSURE 
LOSS 


The  relative  flow  rate  inside  the  flame  tube  of  the  I  cross- 
(Ctlon  combustion  chamber  is 


w  u  —  1  —  VtiAjV,At  —  1  —  v,i<aati 
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From  Eq.  (3-1)  and  (3-3),  we  know  that  we  only  need  to  know 
the  value  In  order  to  compute  the  value.  It  is  not  hard 
to  see  that  what  the  value  of  represents  is  the  average  flow 
rate  coefficient  of  all  the  orifices  behind  the  ith  a^ray  of  orifices. 
The  results  of  calculation  and  experiment  indicate  that  the 
value  may  be  expressed  as  a  function  of  the  X$  value,  where 


<A  — 


/  =»  i 


(3-*D 

(3-5) 


Figure  7  shows  the  cj  curve  obtained  from  the 
experimental  values.  The  horizontal  coordinate  in  the  plot  is 
X<|>  and  the  vertical  coordinate  is  «,  -  ebii  •  w  .  In  the 

Eq.  emax  is  the  coefficient  determined  by  <j>X_Q  ,  see  Figure  8. 

luaX 

Notice  that  the  air-inlet  orifice  in  the  head  section  of  the  flame 


tube  is  not  included  in  the  value  of  Xmax-  From  the  diagram  it 
can  be  seen  that  all  the  experimental  points  of  various  a^^  and 
Xmax  values  basically  fall  on  one  curve.  Thus  curve  may  be 
represented  as  • 


"•  **  °  «^(l  -  ,-».***■> 


(3-6) 


Figure  7a.  Relation  between  ug  and  X  (for  small  X$  values). 
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Figure  7b.  Relation  between  oj£  and  <j>X  (for  large  X4»  values) 


The  emax  curve  in  Figure  8  may  be  approximated  as 


0.725  +  0.090V  -  1 


(3-7) 


From  the  definition  of  Eq.  (3-2),  we  can  see  that:  u  =0,  for 

V  A 

the  last  array  of  orifice  cross-section  of  the  flame  tube,  i.e. , 
the  outlet  cross-section  of  the  combustion  chamber,  and  1  for 

the  cross-section  of  the  first  row  of  orifices  of  the  flame  tube 

when  ♦W"- 

Since  the  coordinate  parameter  of  the  w  curve  is  w£  and 
not  directly  w,  hence  it  is  more  convenient  to  write  (3-1)  in 
the  following  form 


wu  -  l  - 


(3-8) 
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Figure  8.  Relation  between  e  and  <f>X 
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The  a)  curve  in  Figure  7  is  obtained  from  the  data  for  cooling 
air  state  (water  flow  simulation)  when  there  is  no  vortex  generator 
on  the  test  piece.  However,  the.  large  quantity  of  computational  results 
indicate  that  the  w  curve  in  Figure  7  is  truly  a  meaningful  curve 
in  the  general  sense.  It  is  capable  of  taking  into  account  the 
effects  of  many  factors  such  as  J  l  \  A  max  \  .‘i  hO  and  combustion  heat 
release,  etc.,  and  for  flame  tube  or  combustion  chamber  with  variable 
cross-sectional  area,  the  w  curve  above  and  some  of  the  computational 
formula  (including  the  dynamic  equation)  to  be  introduced  below  are 


still  applicable.  Their  degree  of  accuracy  is  even  surprising.  We 
shall  illustrate  this  point  below  by  using  as  a  concrete  example 
a  staged  flame  tube.  The  combustion  chamber-  has  constant  cross- 
section  (fla.  +  au  —  const)  .  The  flame  tube  is  staged  with  6  rows 
of  orifices  (all  with  flat  surface  <?/. ,  —  const  —  0.1333)  )  except  the 

vortex  generator  (<-,  —  0.60,  aM  —  o.os)  .  For  the  two  front 

rows,  a,  —  0.5  ,  the  middle  two  rows,  «/ **  0.6,  and  the 

last  2  rows  a/  — 0.7  .  Heating  ratio  is  0=*l  .  The 

comparison  of  the  calculated  results  by  the  average  flow  rate  coeffi¬ 
cient  method  to  the  results  in  Reference  [31  (see  Reference  [3] 

Figure  lMa  and  Figure  23)  is  shown  in  the  table  below. 


wi 


A 

0 

0.167 

0.333 

0.500 

O.667 

0. 831+ 

1.000 

Average 
flow  rate 
coefficient 
method 

0.116 

0.212 

0.369 

0.1+82 

O.6U5 

0.802 

1.000 

Ref. [3] 

0.125 

0.231 

0.363 

0 . 494 

0. 6i+o 

0.802 

1.000 

The  at  curve  obtained  above  is  based  on  the  experimental  data 
of  the  inlet  orifice  in  this  experiment.  The  relation  between  the 
flow  rate  coefficient  c°  (the  superscript  "o"  represents  the  inlet 
orifice  in  our  experiment  here)  and  the  flow  parameter  Fp  is  shown 
in  Figure  1.  When  the  geometric  shape  of  the  inlet  orifice  we  want 
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to  calculate  Is  different  from  the  inlet  orifice  used  in  this  experi¬ 
ment  so  that  its  flow  rate  coefficient  C  and  F  relation  is  different 

P 

from  the  relation  in  Figure  1,  the  effect  of  the  variation  in  flow 
rate  coefficient  must  be  taker  into  account.  Under  this  circumstance, 
we  may  convert  the  various  shaped  inlet  orifice  area  Ah  to  a£  according  /I 
to  the  principle  of  CAh—  constant,  i.e.,  the  A°  in  Eq.  (3-5)  should 
be  calculated  by  the  equation  below 

Al  -  (c/c°)J,  (3-9) 

Results  of  calculation  show  that  the  flow  rate  coefficient  of  the 

inlet  orifice  is  basically  only  related  to  xV  ;  hence  the  value 

of  the  flow  parameter  Fp  wll  also  be  determined  basically  by  the  value 

of  xVl>  •  The  relational  curves  between  c°  and  Fp  and  Xy/1>  may 

be  formed  from  our  experimental  result  as  shown  in  Figure  9  and  Figure 

o 

10.  According  to  Figure  9,  the  c  value  may  be  determined  and  the 
value  of  F  may  be  determined  from  Figure  10  and  hence  the  c  value  for 
other  inlet  orifice  shapes  may  be  found  from  the  given  C  — 
relationship  (this  relation  must  be  given  first  in  the  calculation 
of  combustion  chamber  flow  rate  distribution). 

As  mentioned  above,  it  is  necessary  to  know  c/c°  when  X  is 
to  be  determined  and  the  value  of  c/c°  is  to  determined  based  on  the 
X  value.  Hence 


Figure  9a.  Relation  between 
small  orifice  flow  rate  coef¬ 
ficient  c°  and  (for 

small  xv'j  values). 
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Figure  9b.  Relation  between 
small  orifice  flow  rate 
coefficient— c°  and  W7- 
(for  large  x^  values).  •<uo 
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Figure  10a.  Relation  between 
flow  parameter  F  and  xS* 
(for  small  xvT?  p  values). 


Figure  10b.  Relation  between 
flow  parameter  F  and 
(for  large  p  values). 


A  computational  process  is  required  in  the  calculation.  Under 
these  circumstances,  it  is  best  to  calculate  backward  from  the 
last  row  of  orifices  of  the  flame  tube  toward  the  front. 


In  the  section  below,  we  shall  discuss  the  method  of  finding 

V.  and  V.  values: 

|  Jo  jex 

The  value  of  the  jet  velocity  v^Q  at  the  initial  flame  tube 

cross-section  may  be  easily  determined  from  Eq.  (3-3).  However, 

the  value  of  _ needs  to  be  found  from  the  curves  based  on  the 

■  max 

!  «  value  of  X<)>  and  the  value  of  $  is  still  unknown.  Hence,  in  the 

calculation,  one  may  first  take  a  $  value  (e.g.  for  cooling  air 
we  may  take  <t>  1  and  when  heating  ratio  >1,  we  may  choose  larger 

values  for  <p  ),  find  u>max  from  w  curve,  then  calculate  VjQ  value 
from  Eq.  (3-3).  We  shall  see  below  that  the  jet  velocity  Vjex  at 
the  outlet  cross-section  of  the  combustion  chamber  may  be  calculated 
from  the  value  of  VjQ,  and  therefore  a  new  value  of  may  be  found. 

Based  on  this  <p  value,  we  can  re-calculate  VjQ  until  the  2  successive 
values  of  VjQ  are  close  to  each  other. 

Table  2  gives  the  comparison  between  the  VjQ  value  as  calculated 

above  and  the  experimental  results.  We  can  see  that  the  VjQ  values 
calculated  from  the  u>  curves  are  very  close  to  the  experimental 
values. 
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Table  2.  Comparison  of  the  results  calculated  from  the  average 

flow  rate  coefficient  method  to  the  experimental  values. 


te8t 

'  method*^P^*ce 
para-- \ 
meter  \. 

a,  -0.79 

v-«- 

6.52  3 

t,  —0.79 

•V.„  - 

5.415 

»,  —0.6 

.H  — 

6.523 

_ 

j,  -0.63 

^.U  — 

5.435 

_ 

d,  —0.63 

v.„- 

4.348 

>,-0.63 

Af.-- 

2.174 

a,  -0.63 
-v.„  - 

1.087  j 

lit  —0.43 

AC  .ii  — 

6.523 

Lt, -0.43 

v.„- 

5.435 

average  flow 

Hi 

mm 

race  coeffi. 
method 

5.25 

5.26 

B 

3.35 

4.74 

2.02 

1.86 

■Hi 

■ 

1 

■B 

| expert,  value 

m 

2.97 

2.91 

na 

|  3.40 

4.63 

1.94 

1.92 

■  n 

ss 

mtm 

mm 

■ 

Efl 

1.03 

0.96 

0.97  , 

S3 

1.29 

jig 

BBS 

1 

■ 

ES 

■ 

B 

expert,  value 

IB 

■ 

1.16 

a 

1.40 

2.62 

1.34 
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It  is  necessary  to  know  the  value  of  4>  in  order  to  make  use 
of  the  a)  curve  in  Figure  7  in  computations.  It  should  be  pointed  out 
that  the  vjex  value  in  the  horizontal  coordinate  $X  of  Figure  7  has 
not  taken  into  account  the  total  pressure  loss  of  the  2  passages. 
Detailed  calculation  shows  that  the  effect  of  whether  or  not  taking 
into  account  of  the  v.  of  the  total  pressure  loss  of  the  2  passages 

v  X 

in  getting  the  w  curve  in  the  calculation  of  flow  rate  distribution 
and  total  pressure  loss  is  basically  the  same,  but  the  latter  is  much 
simpler.  Under  this  condition,  we  have 


P*l  Pie* 


(3-10) 


From  the  flame  tube  dynamic  equation: 
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We  get 
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under  the  condition  of 


—  con*t  and 
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,  then 
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The  combustion  chamber  friction  coefficient  is 


e  _  ga  ~  £as  - 

p.vi/1 


(3-1^) 


It  should  be  pointed  out  that  the  average  flow  rate  coefficient 
method  introduced  above  is  also  applicable  to  the  situation  of  ring 
type  combustion  chamber.  Due  to  the  limitation  of  space,  we  shall  not 
discuss  it  here. 


4.  -COMPARISON  OF  EXPERIMENTAL  RESULTS  AND  COMPUTATIONAL  RESULTS 
FOR  COMBUSTION  CHAMBER  FLOW  RATE  DISTRIBUTION  AND  TOTAL 
PRESSURE  LOSS 


In  Figures  11,  12,  5  and  6  are  presented  respectively  the  com¬ 
parison  of  the  calculated  results  and  experimental  results  for  flow 
rate  distribution  and  total  pressure  loss.  From  the  diagrams,  it 
can  be  seen  that  the  results  calculated  by  solving  the  one  dimensional 
flow  equations  (not  taking  into  account  of  total  pressure  loss  of 
2  passages)  and  the  results  calculated  by  the  average  flow  rate 
coefficient  are  in  good  agreement  with  the  experimental  result  while 
those  calculated  by  the  flow  resistance  method  differ  more  from  the 
experimental  result. 


Figure  11.  Comparison  between  Figure  12.  Comparison  between 

calculated  results  of  flow  the  computed  values  of  flow 

distribution  and  experimental  rate  and  the  experimental 

value.  o»,-<n79,  6.123)  value  (*,^0.43,  xmm-6.riT) 

1  -  Experimental  results;  1  -  Experimental  results; 

2  -  calculated  results;  2  -  calculated  results; 

3  -  average  flow  rate  3  -  average  flow  rate 

coefficient;  4  -  Knight  method  coefficient;  4  -  Knight 

method 


SYMBOL  EXPLANATION 


■ 

A, a  -  represent  respectively  the  x-i/c,  - 


area  and  relative  area  (A/Ar) 

<?  “  2  Ak*I  2  A"i 

i.l  I- 1 

c  -  orifice  flow  rate  coefficient 
(corresponding  to  total  static 
pressure  difference) 

F,~(vt/v.y  _  fiow  parameter 

”  area  parameter 

i  -  # 

n  -  number  of  orifices  (not  including 
the  orifice  in  the  head  of  the 
flame  tube) 

p,pmy^p*  -  respectively  static 

pressure,  total  pressure 
and  total  pressure 
difference 

-  weight  flow  rate  through 
flame  tube  orifice 

i  -  (0/2)1'’  -  velocity  head 

V,v  -  respectively  velocity  and 
relative  velocity  (V/Vr) 

I 

-  2  AQil*iA'v' 

t  m  I 

aw  -  ao/o,  -  relative  flow  rate 
through  flame  tube 
orifice 


5-(»r.  y‘  ~  fiction 

coefficient 

SUBSCRIPT 

a,l  -  represent  respectively 

2  passage  and  flame  tube. 

h  -  flame  tube  orifice 

i  -  order  number  of  the  cross- 
section  of  the  orifice 
array  or  in  front  of  the 
orifice  array  of  the  flame 
tube 

j  -  order  number  of  the  orifice 
jet  or  orifice  array 

o  -  head  orifice  of  the  com¬ 
bustion  chamber  (including 
the  vortex  generator) 

ex,  in  -  respectively  the 

outlet  and  inlet  cross- 
section  of  the  combus¬ 
tion  chamber  or  the  2 
passage  way. 

r  -  reference  cross-section 
(taken  as  a,  -  a.  +  a,  in 
this  paper) 


SUPERSCRIPT 

o  -  corresponds  to  the  inlet 
orifice  used  in  our  experi¬ 
ment 

*  -  stagnation  parameter 
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